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Introduction 


Conventional  mammography  is  limited  in  its  sensitivity  and  specificity  for  breast  cancer 
detection,  since  it  relies  on  the  small  difference  in  x-ray  attenuation  between  breast  tissues  and 
the  lesions.  As  x-ray  wave  passes  through  breast,  x-ray  undergoes  phase-shifts  as  well.  X-ray 
phase-shift  differences  between  tissue  and  lesions  are  several  hundreds  times  larger  than  their 
attenuation  differences.  Hence  the  phase-sensitive  breast  imaging  has  the  potential  to  greatly 
enhance  the  cancer  detection  sensitivity  and  specificity  and  reduce  the  radiation  doses  to 
breast.  It  is  challenging  to  implement  phase-sensitive  x-ray  breast  imaging  because  one  should 
provide  the  partially  coherent  x-ray  illumination,  maintain  short  exposure  times,  and  limit 
radiation  doses.  Moreover,  in  order  to  fully  exhibit  tissue  phase  contrast,  one  has  to  perform 
phase  retrievals.  But  the  most  commonly  used  phase  retrieval  method  is  unstable  to  image 
noise  and  requires  high  radiation  dose  levels  to  work,  thus  this  common  phase  retrieval  method 
is  not  feasible  in  clinical  imaging  where  radiation  dose  control  is  critically  important. 

Taking  up  these  challenges,  the  long-term  objective  of  the  project  is  to  develop  a  low- 
dose  and  quantitative  x-ray  phase  imaging  technique  for  facilitating  breast  cancer  detection. 

The  Specific  Aims  of  the  project  are:  (1 )  Develop  a  prototype  phase-imaging  system  enabling 
the  phase  retrieval,  that  is,  the  reconstruction  of  objects  phase-maps.  The  system  hardware 
comprises  a  micro-focus  tube  operating  at  high  tube  voltages,  a  high-resolution  photostimulable 
phosphor  based  detector  (CR-system).  An  innovative  phase  retrieval  method  will  be  developed 
to  retrieve  tissue  phase  map  from  a  single  recorded  image.  (2)  Validate  the  accuracy  of  the 
reconstructed  tissue  projected  electron  densities;  validate  the  many-fold  radiation  dose 
reduction  to  be  achieved  with  the  proposed  system;  conduct  subjective  measurements  to 
characterize  the  performance  of  the  proposed  system. 


Body 


Toward  to  the  Specific  Aims  set  as  above,  we  have  accomplished  following  works 
associated  with  each  task  outlined  in  the  approved  Statement  of  Work. 

Development  of  a  prototype  x-ray  phase  imaging  system  (Task  1) 

Design  and  build  the  system  prototype 

In  this  project  we  designed  and  built  a  unique  x-ray  imaging  prototype  using  high  tube- 
voltages  for  implementing  the  phase  contrast  imaging  and  enabling  radiation  dose  reduction. 
The  key  design  issue  for  an  inline  phase  contrast  imaging  system  is  how  to  optimize  the  phase 


Figure  1.  Relative  phase  contrast 

visibility  |  RPF(u)\  varies  for  different  Ri  and 

R2  of  a  system  with  a  SID  =1 .5  m.  In  this 
plot  we  assumed  a  uniform  circular  focal 
spot  of  50  pm  in  diameter,  and  a  detector 
pitch  of  43.75  pm ,  60.5  keV  x-ray  and  a 
fine  object  structure  of  15  line-pairs/mm. 
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contrast  visibility  by  balancing  the  conflict  requirements  on  x-ray  spatial  coherence  and 
reasonable  exposure  times  .  In  our  design  we  tackled  this  issue  by  using  a  method  based  the 
relative  phase-visibility  factor  (RPF),  a  phase  imaging  metric  that  we  previously  developed  [Wu, 
Liu  2004].  We  computed  the  RFP  as  a  function  of  different  configurations  with  varying  x-ray 
wavelengths,  tube  focal  spot  sizes,  different  source-object-  and  object-detector-distances 
(denoted  by  R ^  and  R2  respectively),  finest  object-structure  sizes,  and  the  detector  sampling 
pitches.  Form  the  computed  RPFs  we  determined  the  optimal  system  configuration  for  phase 
contrast  visibility.  Figure  1  shows  the  RFP  as  a  function  of  Ri  and  R2  for  a  system 
configuration.  Figure  1  shows  that  the  optimal  magnification  factors  are  found  to  be  in  the  range 
of  2  to  2.3.  The  experiment  with  various  phantoms  (such  the  step-edge  phantom, 
mammography  BR3D  phantom,  contrast  details  phantom,  etc.)  confirmed  optimal  phase 
contrast  visibility  with  the  designed  system  configurations. 

According  to  these  optimal  system  configurations 
determined  from  above  RPF-analysis,  our 
subcontractor  Dr.  Liu  and  his  group  integrated  the 
hardware  components  successfully  built  a 
prototype  system  as  is  shown  in  Fig.  2.  In  order  to 
operate  the  prototype  in  both  conventional  and 
phase  contrast  imaging  modes  for  comparisons, 
the  imaging  and  measurement  components  (the  x- 
ray  tube,  CR-plate)  are  mounted  on  an  optical  rail, 
which  allows  to  vary.  Note  that  the  object  -detector 
distance  R2  is  kept  relatively  large  as  determined 
by  the  RFP-analysis  to  allow  phase-shifted  x-rays 
diffract  to  form  phase-contrast  images.  The  source 
is  a  micro-focus  tube  (Model  L8 121-01, 

Hamamatsu  Photonics)  with  integrated  high 
voltage  power  supply.  The  tube  has  a  tungsten- 
target  and  a  uniform  circular  focal  spot  of  50  |im  in 
diameter.  The  mounted  detector  is  a  Konica  Minolta  Konica  large  FOV  (24  x  30cm) 
mammography  photostimulable  phosphor  plate  and  a  Konica  REGIUS  1 90  Reader  of  a 
43.75  nm  sampling  pitch.  The  system  has  been  used  in  acquiring  phase  contrast  projection 
images  of  various  phantoms  in  this  project. 


Fig.  2.  A  high-energy  x-ray  phase  contrast 
imaging  system 


Develop  the  phase  retrieval  methods 

With  a  micro-focus  x-ray  tube  the  achievable  x-ray  spatial  coherence  is  limited,  hence  a 
x-ray  tube-based  system  can  generate  only  moderate  phase-contrast  at  tissue  interfaces  and 
boundaries  with  dark-bright  diffraction  fringes.  In  order  to  fully  exhibit  the  tissue  phase  contrast, 
one  should  perform  the  phase  retrieval  from  the  acquired  phase-contrast  images.  The  phase 
retrieval  is  a  procedure  of  retrieving  the  phase-shift  map^(r  )  of  an  object  from  its  phase- 

sensitive  projections,  where  cp(r)  =  -Xre  J  pe(r ,s)ds  and  pe (r, 5)  denotes  the  tissue  electron 

ray 

density,/;  =  2.818  x  1CT15  m  is  the  classical  electron  radius  and  1  is  x-ray  wavelength.  By 

means  of  phase  retrieval  from  phase-contrast  images,  one  can  disentangle  the  phase  contrast 
from  mixed  attenuation-  and  phase-contrast  in  phase-contrast  images.  Moreover,  a  retrieved 
tissue  phase  map  provides  a  map  of  imaged  object’s  projected  electron  densities  (the  ray 
integrals  of  electron  densities)  as  well.  The  challenge  for  developing  breast  x-ray  phase  imaging 
is  how  to  overcome  the  difficulty  for  implementing  phase  retrieval  with  breast  imaging. 
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In  general  one  needs  to  acquire  at  least  two  and  mostly  four  to  five  images  with  different 
object-detector  distances.  This  requirement  is  imposed  by  the  x-ray  wave  nature  and  its  wave 
equation.  The  most  commonly  used  method  in  the  literature  for  phase  retrieval  in  the  literature  is 
based  on  the  Transport  of  Intensity  Equation  (TIE)  [Allen,  Oxley.  2001],  According  to  the  TIE- 
based  method,  one  needs  to  acquire  at  least  two  images,  one  contact  radiograph,  and  one 
phase  contrast  image,  for  retrieving  a  phase  map  of  an  object.  While  the  TIE-based  phase 
retrieval  method  is  computational  effective  for  the  cases  with  strong  phase  contrast  effects  and 
low  noise  levels,  we  found  that  it  will  fail  for  phase  retrieval  when  the  acquired  projection  images 
are  contaminated  with  relatively  strong  noises.  In  other  words,  we  found  that  the  common  TIE- 
based  method  is  not  robust  against  noise.  Suppressing  x-ray  quantum  noise  may  require  using 
high  radiation  doses  in  the  image  acquisitions.  For  clinical  applications  it  is  imperative  to  limit 
and  reduce  radiation  doses  involved,  hence  it  is  critical  to  develop  low-dose  robust  phase 
retrieval  methods  for  clinical  applications  of  phase  contrast  imaging,  as  in  clinical  applications 
where  radiation  dose  constraints  dictates  relatively  high  noise  levels.  Moreover,  the  TIE-based 
method  needs  to  acquire  multiple  (at  least  two)  images  with  different  object-detector  distances. 
This  multiple-exposure  requirement  adds  to  the  radiation  dose  burden  and  makes  the  procedure 
cumbersome.  To  overcome  these  difficulties  in  implementing  phase  retrieval  in  clinical  imaging, 
we  developed  two  phase  retrieval  methods  as  described  in  following: 

(1)  Development  of  the  Phase-Attenuation  Duality  based  (PAD-based)  phase 

retrieval  method 

First,  we  found  that  when  high  energy  x-rays  are  employed,  we  can  robustly  retrieve  the 
map  of  a  soft  tissue  object  from  just  a  single  phase-contrast  projection  of  the  object.  In  fact, 
analyzing  the  tissue  attenuation  data,  we  realized  that  the  soft  tissues’  attenuation  cross 
sections  can  be  approximated  by  that  of  x-ray  incoherent  scattering,  as  long  as  x-ray  energies 
are  of  about  60  keV  or  higher.  Under  this  circumstance,  we  found  that  soft-tissue  attenuation 
and  soft-tissue  phase-shifts  are  all  related  to  the  projected  electron  density.  We  call  this 
correlation  as  the  phase-attenuation  duality.  Utilizing  this  notion  of  the  phase-attenuation  duality 
(PAD),  we  were  able  to  simplify  x-ray  wave  equation,  and  discovered  an  innovative  method  to 
retrieve  a  phase-map  from  just  a  single  recorded  image  [Wu,  Liu,  Yan  2005],  We  call  this  novel 
method  as  the  PAD-based  phase  retrieval  method.  With  this  method  the  cumbersome 
procedures  associated  with  required  exposures  for  the  phase  retrieval  are  eliminated,  and  more 
importantly,  the  radiation  doses  involved  are  reduced  as  well. 

In  the  project  we  employed  high  tube  voltages  such  as  120  kVp-150  kVp  for  breast 
phantom  imaging  for  shortening  the  exposure  times  and  reducing  radiation  doses.  Phase 
retrieval  from  a  single  phase-contrast  image  forms  a  key  strategy  for  this  project.  However,  we 
noted  that  our  previous  PAD-based  phase  retrieval  method  becomes  inaccurate  in  the  cases  of 
high-resolution  and  large  object-detector  distances.  Extending  our  previous  results,  in  this 
project  we  developed  a  more  accurate  PAD-based  phase  retrieval  method.  We  noticed  that  the 
Fresnel  diffraction  can  be  thought  as  a  x-ray  wave-front  filtering  with  the  Fresnel  propagator 
exp(/;zvtzM2),  where  A  is  x-ray  wavelength  and  u  is  the  wave-front's  transverse  spatial 
frequency  vector,  and  z  is  the  object-detector-distance.  So  the  phase  of  the  propagator,  i.e. 
nXzif ,  determines  how  strongly  the  exiting  wavefront  will  be  diffracted  in  the  propagation.  The 
previous  PAD-based  phase  retrieval  formula  was  derived  by  linearizing  Fresnel  propagator 
phases  [Wu,  Liu,  Yan  2005],  Therefore,  in  order  to  develop  robust  phase  retrieval  algorithms 
with  high  accuracies,  one  should  not  linearize  Fresnel  propagators,  rather  one  should  consider 
the  full  exponential  propagators  instead.  In  this  project  we  derived  a  new  PAD-based  phase 
retrieval  method  using  the  full  exponential  propagators.  This  new  phase  retrieval  method  lays  a 
sound  foundation  for  implementing  high-resolution  phase-sensitive  imaging  of  soft  tissue 
objects.  For  details,  please  refer  to  [Wu  and  Yan  2009]. 
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(2)  Development  of  the  Attenuation-Partition  based  (AP- 
based)  phase  retrieval  method 

The  PAD-based  phase  retrieval  method  described  above 
can  be  used  only  for  soft  tissue  body  parts  imaged  with  high- 
energy  x-rays.  One  still  needs  to  develop  a  robust  phase 
retrieval  method  for  general  clinical  imaging  applications  as  well. 
Extending  the  PAD-method  to  general  cases,  we  developed  the 
attenuation-partition  based  (AP-based)  method  for  general 
phase  retrievals.  In  this  method,  we  utilize  the  correlation 
between  the  x-ray  phase  shift  and  its  attenuation  to  eliminate  the 
intrinsic  singularities  involved  in  the  TIE-based  phase  retrieval 
method.  These  singularities  may  amplify  the  noise  in  the 
acquired  images  and  result  in  instability  in  phase  retrievals.  In 
the  AP-based  method,  we  first  extract  the  Compton  scattering 
cross-sections  from  the  attenuation  image  and  retrieve  an 
approximate  phase  map.  We  then  iteratively  incorporate  the  rest 
of  the  attenuation  cross-section  into  phase  retrieval  by 
repeatedly  correcting  errors,  which  is  calculated  as  the 
differences  between  the  measured  phase-contrast  image 
intensities  and  the  Fresnel  diffraction  estimates  [Yan,  Wu,  Liu 
2010],  We  first  proposed  the  idea  in  2008,  and  in  this  project  we 
modified  the  original  iterative  algorithm  coding  for  speeding  the 
iteration  convergence,  and  performed  quantitative  and 
systematic  evaluation  of  its  robustness  and  efficiency.  We  systematically  compared  the 
performance  of  this  new  algorithm  with  two  widely  used  phase  retrieval  algorithms  in  the 
literature  including  the  TIE-based  method  described  above.  We  showed  that  the  AP-based 
method  converges  faster  than  the  Gerchberg-Saxton  algorithm  in  the  Fresnel  diffraction  regime, 
and  is  much  more  robust  against  image  noise  than  the  TIE-based  method.  For  details,  please 

Moreover,  we  also  analyzed  the  phase  maps  retrieved 
from  the  experimental  images  of  a  phantom  by  using  the  AP- 
based  and  TIE-based  methods,  respectively.  We  applied  these 
two  methods  to  experimental  projection  images  of  an  air-bubble 
wrap  phantom  for  retrieving  the  phase  map  of  the  bubble  wrap. 
Fig.  3-(upper)  shows  the  retrieved  sample's  phase  map  using 
TIE-based  method.  Apparently  the  retrieval  with  the  TIE-based 
method  is  very  unsatisfactory,  as  the  bubble  rims  are  hardly 
visible  in  the  extremely  cluttered  background  in  the  phase  map. 
In  a  stark  contrast,  Fig.  4-(upper),  which  was  retrieved  using  the 
AP-method,  depicts  the  bubble  rims  prominently.  In  order  to 
gauge  the  accuracies  of  the  retrieved  phase  values,  we 
measured  the  thicknesses  of  the  flat  bases  of  the  wrap,  and 
found  that  the  approximate  phase-shifts  at  the  flat  bases  of  the 
wrap  is  about  -4.2  radian.  On  the  other  hand,  according  to  the 
phase  profile  in  with  the  AP-based  method  (Fig.  4(lower)),  the 
average  phase  shift  at  the  flat  bases  of  the  wrap  is  -1 .05  radian, 
the  two  phase-shift  estimates  are  reasonably  close.  In 
comparison,  the  phase  profile  with  the  TIE-based  method  (Fig. 
3(lower),)  depicts  just  messy  up-down  peaks  buried  in  cluttered 
and  noisy  background.  According  to  this  profile,  the  phase-shift 
values  at  the  bases  fluctuate  over  a  wide  range  from  0  to  +87.5 


refer  to  [Yan,  Wu,  Liu  2010], 


Fig.  4.  (Upper)  Retrieved  phase 
maps  of  the  bubble  wrap  with 
the  AP-based  method.  (Lower) 
Profile  of  the  retrieved  phase 
values  along  the  marked  line. 


Fig.  3.  (Upper)  Retrieved  phase 
maps  of  the  bubble  wrap  with 
the  TIE-based  method  and 
mended  with  the  Tikhonov 
regularization.  (Lower)  Profile 
of  the  retrieved  phase  values 
alono  the  marked  line. 
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radian  with  an  average  of  +27.3  radian,  contradicting  to  the  fact  that  x-ray  phase  shifts  should 
be  negative.  This  stark  comparison  shows  that  the  conventional  TIE-based  phase  retrieval 
method  is  unstable  against  the  noise  in  the  wrap’s  projection  images,  while  the  AP-based  phase 
retrieval  method  is  shown  in  the  experiment  to  be  superior  to  the  TIE-based  method  for  the 
robustness  in  performing  the  phase  retrieval.  For  details,  please  refer  to  [Yan,  Wu,  Liu  2011], 

Extending  the  research  works  on  the  phase  retrieval  methods,  as  a  by-product  of  this 
project,  we  made  a  progress  in  understanding  the  apparent  linear  attenuation  coefficients  in  the 
inline  phase  contrast  tomography.  In  the  inline  phase  contrast  x-ray  tomography  the 
reconstructed  apparent  linear  attenuation  coefficient  values  in  the  tomography  may  be  greatly 
larger  than  sample's  linear  attenuation  coefficients  or  even  be  negative.  These  "artifacts"  may 
cause  faulty  interpretation  of  sample  structures,  and  impede  even  qualitative  characterizations 
for  tissues  and  materials.  In  this  period  we  derived  a  general  formula  to  quantitatively  relate  the 
apparent  linear  attenuation  coefficient  values  in  cone-beam  phase  contrast  tomography  to 
sample's  linear  attenuation  coefficients  and  refractive  indices.  This  formula  overcomes  the 
gross  inaccuracy  of  the  existing  formula  in  the  literature  in  analyzing  high-resolution  phase 
contrast  tomography,  and  it  will  be  useful  for  correctly  interpreting  and  quantifying  the  apparent 
linear  attenuation  coefficients  in  cone-beam  x-ray  phase  contrast  tomography  in  the  biomedical 
and  material  science  applications.  This  result  advanced  the  understanding  of  the  image  contrast 
mechanism  of  phase  contrast  tomography.  For  details,  please  refer  to  [Yan  and  Wu  2011], 

Development  of  system  software  package 

In  this  project  we  developed  efficient  coding  with  the  software  modules  for  the  image 
registration,  flat-filed  correction,  and  phase  retrievals.  In  addition,  we  developed  as  well  phase 
contrast  imaging  simulation  software  incorporating  parallel  programing. 


Characterization  of  the  x-ray  phase  imaging  system  (Task  2) 

Investigating  the  effects  of  x-ray  beam  hardening  on  detective  quantum  efficiency 
(DQE)  and  radiation  dose 

A  comprehensive  quantitative  determination  of  the  performance  of  an  imaging  system  is 
the  detective  quantum  efficiency  (DQE),  as  it  represents  the  transfer  function  of  the  signal-noise 
ratio  of  the  system.  The  goal  of  this  study  was  therefore  to  provide  an  initial  investigation  of  the 
effects  of  a  range  of  beam  hardening  levels  on  the  DQE  and  the  average  glandular  dose.  The 
DQE  and  the  average  glandular  dose  were  both  measured  under  the  same  experimental 
conditions  for  a  range  of  beam  hardening  levels.  The  comparison  of  all  quantities  comprising  the 
DQE  exhibits  very  close  correlation  between  the  results  obtained  without  added  beam 
hardening  to  the  results  corresponding  to  the  range  of  beam  hardening  levels.  For  the  specific 
experimental  conditions  utilized  in  this  study,  the  results  are  an  indication  that  the  use  of 
hardened  x-ray  beam  holds  the  potential  to  reduce  the  radiation  dose  without  decreasing  the 
performance  of  the  system.  For  details,  please  refer  to  [Wong,  Wu,  Liu.  2011], 

Evaluating  Image  quality,  phase  map  accuracy  and  associated  radiation  dose 
reduction 

We  have  conducted  preliminary  imaging  experiments  using  100,  120,  140  kVp  and  with 
four  different  phantoms,  the  ACR  MAP  phantom,  a  contrast-detail  phantom,  an  acrylic  step- 
edge,  and  a  breast  tissue-equivalent  phantom.  As  current  breast  imaging  (mammography  and 
breast  tomosynthesis)  employs  exclusively  low  x-ray  tube  voltages  (25-40  kVp),  our  strategy  is 
to  utilize  the  phase  contrast  extracted  out  with  the  phase  retrievals  to  compensate  for  the 
attenuation-contrast  loss  associated  with  the  use  of  high  tube-voltages,  thereby  achieving 
significant  radiation  doses  reduction.  We  images  the  phantoms  with  different  geometric 
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configurations 
including  the 
settings  with 
magnification 
factors  of  1 ,  2,  2.5, 
3.  Observer 
studies  were 
performed  utilizing 
the  ACR  MAP  and 
contrast-detail 
phantoms  for 
qualitative  and 
quantitative 
comparisons,  such 
as  those  in  terms  of 

Fig.  5.  Two  images  of  an  ACR  accreditation  phantom  acquired  with  the  100  kVp  and  ^he  feature  scoring 

the  same  dose  to  phantom,  (left)  Phase  contrast  image  of  the  phantom  with  M  =  2.  and  contrast-noise 

(right)  Contact  radiograph  of  the  phantom.  ratio  (CNR).  In  this 

way  we  validated 

the  optimal  geometric  configurations  for  phase  contrast  visibility  as  determined  in  the  system 
design  based  on  the  theoretical  phase  contrast  visibility  analysis  [Wu,  Liu  2004],  In  addition  to 
phantom  image  quality  evaluation,  in  order  to  characterize  the  performance  of  the  prototype 
system,  we  checked  as  well  the  accuracy  of  the  retrieved  phantom  phase  maps,  and  verified  the 
radiation  dos  reduction  achieved. 

Study  on  ACR  MAP  Phantom 


Fig.  6.  Phase  map  of  an  ACR  accreditation 
phantom  retrieved  by  using  the  duality-based 
method  from  a  single  phase-sensitive  projection. 
The  projection  was  acquired  with  120  kVp  and  a 
glandular  dose  of  85  mrad,  a  very  low  dose 
compared  to  that  employed  in  mammography.  The 
phase  map  provides  as  well  the  phase-shift 
values,  e.g.,  the  phase-shift  of  the  first  mass 
including  the  acrylic  base  is  found  as  -1106 
radians. 


Fig.  5  shows  two  images  of  a  ACR  MAP 
phantom  acquired  with  the  prototype  system. 

The  Fig.  5-(left)  is  a  phase-sensitive  projection 
acquired  at  100  kVp  and  a  magnification  factor  of 
2.  For  comparison,  the  Fig.  5-(  right)  shows  the 
contact  radiograph  of  the  phantom  acquired  with 
the  same  kVp  and  the  same  phantom  dose.  Both 
images  were  acquired  without  any  anti-scatter 
grid.  The  loss  of  the  attenuation  contrast 
associated  with  the  high  tube  voltage  is 
responsible  for  the  poor  contrast  shown  in  Fig.  5- 
(right),  though  the  x-ray  scatter  contributed  to  the 
image  contrast  degradation  as  well.  In  contrast, 
the  phase-sensitive  projection  image  in  Fig.  5- 
(left)  scores  much  better  in  visibilities  of  the 
fibers,  specks  and  masses  in  the  phantom. 

In  order  to  fully  exhibit  the  phase  contrast 
of  phantom  components  and  quantify  their 
projective  electron  densities,  we  retrieved  the 
phantom's  phase  maps.  As  an  example,  shown 
in  Fig.  6  is  a  phase  map  of  an  ACR  phantom 
retrieved  by  applying  our  PAD-based  phase 
retrieval  method  to  a  single  phase-sensitive 
projection  acquired  with  the  prototype  system. 
The  exposure  was  at  120  kVp  and  with  an 
estimated  average  glandular  dose  of  0.85 
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mGy.  As  is  shown  in  Fig.  6,  the  ACR-phantom  phase  map  scores  4  fibers,  3  specks  group,  and 
3.5  masses  for  the  feature-visibilities. 

To  accurately  estimate  the  average  glandular  doses  involved  in  experiments,  we 
extended  our  previously  published  tables  of  the  normalized  average  glandular  doses  (d  N )  for 

conventional  mammography  to  high-energy  x-ray  imaging.  We  performed  the  Monte  Carlo 
simulations  of  x-ray  transport  and  energy  deposition  in  breast  for  photon  energies  ranging  up  to 
150  keV.  The  D„v  -values  for  each  of  photon  energies  were  then  weighted  according  to  the 

measured  x-ray  spectrum  (at  100,  120  and  140  kVp  respectively)  to  computed  theD„v  -values 

associated  with  each  of  the  x-ray  spectra.  In  this  way  we  found  the  average  glandular  dose  is 
the  0.85  mGy  for  the  image  shown  in  Fig.  6.  While  this  dose  is  likely  overestimated  because  of 
the  tube-voltage  ramping  effects,  this  overestimated  dose  value  is  still  very  low  compared  to  the 
phantom  dose  levels  of  1.5-2. 5  mGy  commonly  employed  in  conventional  mammography.  The 
Konica  mammography  CR-plate  detector  employed  in  the  experiment  is  really  not  quantum- 
efficient  for  this  experiment,  because  its  thin  scintillator  layer  is  designed  for  mammography,  not 
suitable  for  high  tube-voltage  imaging  with  such  as  120  kVp  used  in  this  experiment.  More 
contrast-to-noise  ratio  (CNR)  enhancement  and  dose-savings  can  be  achieved  in  future  by 
using  a  high-resolution  detector  with  improved  quantum-efficiencies.  In  fact,  a  high-resolution 
flat  panel  detector  with  higher  quantum-efficiencies  has  been  planned  for  future  experiments. 

Moreover,  the  ACR-phantom  phase  map  in  Fig.  6  provides  as  well  phase-shift  values.  From 
Fig.  6  we  found  that  the  retrieved  phase-shift  of  the  first  mass  including  the  acrylic  base  is 
-1106  radians.  Note  that  x-ray  phase-shift  should  be  of  a  negative  value.  A  more  negative  value 
indicates  a  larger  projected  electron  density.  Compared  to  the  first  mass,  the  retrieved  phase- 
shift  of  the  phantom's  acrylic  frame  differs  by  -35  radians,  and  the  phase-shift  of  the  wax-insert 
(including  the  acrylic  substrate)  differs  by  +10  radians.  To  check  the  accuracies  of  these 
retrieved  phase  values,  we  computed  phase-shifts  based  on  the  ray-tracing  through  a  digital 
model  of  the  ACR  MAP  phantom.  Compared  to  the  ray-tracing  based  phase-shift  value  of  the 
first  mass,  the  phase  value  of  the  acrylic  frame  differs  by  -36  radians,  and  the  phase  value  of 
the  wax-insert  differs  by  +8  radians.  These  results  are  in  excellent  agreement  with  the  phase 
retrieval  results.  However,  the  first  mass's  phase  value  computed  based  on  the  ray-tracing  is 
-887  radians,  which  differs  from  its  retrieved  phase  value  by  20%.  This  difference  may  be  due  to 
the  uncertainty  in  the  average  photon  energy  caused  by  the  tube-voltage  ramping  in  the 
experiment,  and  the  accuracy  can  be  improved  with  better  exposure  control.  In  a  summary,  the 
study  shows  that  the  phase-retrieval  based  high-kVp  phase-sensitive  projection  technique  is 
able  to  significantly  reduce  radiation  doses  and  to  enable  quantitative  imaging  based  on  the 
retrieved  phase-shit  values. 

Extending  this  study,  we  also  found  a  novel  way  to  determine  the  breast  density  by  using 
the  retrieved  phase-shift  values  of  a  breast.  The  volumetric  breast  density  (VBD)  is  the  fraction 
of  the  volume  of  fibroglandular  tissue  in  the  breast.  High  breast  density  inferred  from 
mammograms  is  associated  with  two-  to  four-fold  increase  in  the  risk  for  breast  cancer.  Public 
awareness  of  breast  density  being  as  a  high  risk  factor  pushes  many  states  passing  breast 
density  notification  legislations  for  breast  imaging.  Current  mammography-based  method  for 
VBD-measurement  is  based  on  small  differences  in  the  linear  attenuation  coefficients  of  breast 
fibroglandular  tissue  and  adipose  tissue.  But  breast  linear  attenuation  coefficients  vary  in  a 
complicated  way  with  different  tube  targets,  filtrations  and  operating  voltages.  Hence  one  has  to 
measure  VBD  with  simultaneous  exposure  of  a  calibration  phantom  alongside  the  breast  for 
system  calibrations.  Other  methods  require  laborious  system  calibrations  with  different  tube 
targets/filters  and  voltages.  Note  that  current  digital  breast  tomosynthesis  is  unable  to  provide 
the  VBD  either,  since  the  limited  angular  projections  in  tomosynthesis  cannot  recover  the  breast 
linear  attenuation  coefficients  from  the  image  reconstruction.  Instead  of  measuring  tissues' 
linear  attenuation  coefficients  as  a  surrogate  of  the  breast  density,  we  recently  developed  for  the 
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first  time  an  x-ray  phase-shifts  based  method  of  breast  density  measurement.  In  order  to 
measure  VBD,  we  developed  an  x-ray  phase-shifts-based  method  of  breast  density 
measurement.  We  found  that  the  retrieved  breast  phase  map  from  a  single  projection  can  be 

used  to  compute  the  VBD  asVBD  =  {Zr(-p(r)/(,l)re  -  P^)} /[n pTc(PeJg  - pead)]  , 

where  <p(r)  denotes  the  retrieved  breast  phase  map,  and  Tc\  s  the  breast  compression  thickness, 

p  ,  =  3.448  xlO23/ cm3  and  p  =  3.108  x  1023  /cm3 are  the  electron  densities  of  the  fibroglandular 

tissue  and  adipose  tissue,  respectively.  We  have  numerically  validated  this  novel  method  of 
breast  density  measurement  with  very  good  accuracies.  For  details,  refer  to  [Wu,  Yan,  Liu. 
2012]. 

Study  on  contrast-to-noise  ratios  with  a  contrast  detail  phantom 

We  investigated  the  effect  of  phase  retrievals  on  the  contrast-to-noise  ratios  (CNRs)  with  a 
contrast-detail  phantom  (CD-Phantom).  Shown  in  Fig.  7-(left)  is  the  phase  contrast  projection 
image  of  the  CD-phantom  acquired  at  120  kVp  and  with  a  magnification  factor  of  2.5. 

In  Fig.7-(left)  the  phase  contrast  is  exhibited  as  the  enhanced  edges  of  the  test-disks.  Using  the 
PAD-based  method  as  described  earlier,  we  performed  phase  retrieval  from  the  phantom's 

phase  contrast 
image.  The 
retrieved  phase 
map  of  the  CD- 
phantom  is  shown 
in  Fig.  7-(right).  We 
measured  the  bulk 
CNRs  with  the  disk 
boundaries 
excluded  in  both  the 
images.  We  found 
that  the  disk  CNRs 
in  the  phase  map 
(Fig.  7-(right))  are 
greatly  enhanced, 
and  are  about  four- 
times  the  disk 
CNRs  in  the  phase 

contrast  image  (Fig.  7-(left).)  The  CNR  of  a  given  disk  is  proportional  to  the  square  root  of  the 
radiation  dose  used  in  exposure.  Hence  the  PAD-based  phase  retrieval  has  good  potential  of 
greatly  reducing  radiation  doses  for  achieving  a  given  CNR-level. 

Study  on  the  PAD-phase  retrieval  performance  at  a  synchrotron  facility 

In  a  collaborative  research  with  Dr.  T.  Xiao  group  of  Shanghai  Institute  of  Applied 
Physics,  we  investigated  the  performance  of  the  PAD-based  phase  retrieval  method  in  the 
phase  tomography.  Shown  in  Fig.  8(a-b)  are  the  tomograms  of  a  slice  of  a  four-component 
phantom.  The  image  in  Fig.  8(a)  is  reconstructed  directly  from  900  phase-contrast  angular 
projections,  which  were  acquired  with  60-keV  of  synchrotron  radiation  by  Xioa  group,  using  the 
standard  filtered  backprojection  (FBP)  algorithm.  The  strong  edge-enhancement  demonstrated  in 
Fig.  8(a)  is  resulted  from  the  phase  contrast.  However,  high  noise-levels  are  found  in  this  image  as  well. 
On  the  other  hand,  the  image  in  Fig.  8(b)  is  the  tomogram  that  is  FBP-reconstructed  from  900  angular 
phase-maps,  which  were  retrieved  from  the  same  projection  set  by  using  the  PAD-based  phase  retrieval 
method.  Apparently  the  tomogram  Fig.  8(b)  reconstructed  from  the  angular  phase-maps  demonstrates 


Fig.  7.  (left)  Phase  contrast  projection  image  of  a  contrast  detail  phantom  acquired 
at  120  kVp.  (right)  Retrieved  phase  map  from  the  acquired  phase  contrast  image. 
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Fig .  8.  Reconstructed  tomograms  of  a  sample  slice  at  60  keV.  The  white 
dash-dot  circles  in  the  two  images  are  the  markers,  (a)  Tomogram 
reconstructed  directly  from  900  phase-contrast  angular  projections. 

(b)  Tomogram  reconstructed  from  900  angular  phase-maps  retrieved  by 
using  the  PAD-based  phase  retrieval  method. 


much  higer  bulk  CNRs  for  sample 
components.  Actually  the  bulk 
CNRs  in  Fig.  8(b)  are  about  ten  to 
fifteen  times  higher  than  that  in 
Fig.  8(a).  The  striking  bulk  CNR 
enhancement  achieved  by  using 
the  PAD-based  method  in  this 
experiment  shows  again  that  the 
PAD-based  phase  retrieval 
method  has  good  potential  to 
tremendously  reduce  the  radiation 
doses  involved  in  imaging.  For 
details,  please  refer  to  [Liu  et  al. 
2012]. 


Key  Research  Accomplishments 

•  Designed  and  built  a  high-voltage  x-ray  imaging  prototype  that  enables  x-ray  phase 
imaging  and  radiation  dose  reduction,  as  is  demonstrated  in  phantom  imaging. 

•  Developed  an  innovative  phase-attenuation  duality  based  (PAD-based)  method,  which 
greatly  facilitating  the  low  doe  phase  retrieval  for  soft  tissue  body  parts  such  as  breast. 

•  With  phantom  imaging  we  demonstrated  that  the  PAD-based  phase  retrieval  greatly 
enhance  the  contrast-to  noise  ratios  of  images,  hence  it  has  good  potential  for  reducing 
radiation  doses  in  breast  imaging. 

•  Developed  a  robust  general  phase  retrieval  method,  the  attenuation-partition  based  (AP- 
based)  method  for  general  clinical  imaging.  This  method  was  demonstrated  to  be 
superior  to  the  mostly  used  phase  retrieval  method  in  literature  (TIE-based  method)  in 
performing  robust  phase  retrieval  for  low-dose  imaging. 

•  Advancing  the  understanding  of  the  image  contrast  mechanism  in  phase  contrast 
tomography,  we  derived  a  general  formula  relating  the  apparent  linear  attenuation 
coefficients  to  the  sample's  linear  attenuation  coefficients  and  refractive  indices  in  phase 
contrast  tomography. 

•  Considering  the  high  breast  density  is  one  of  the  biggest  risk  factors  for  breast  cancer, 
we  proposed  for  the  first  time  an  x-ray  phase-shifts  based  method  of  breast  density 
measurement,  as  a  quantitative  imaging  application  of  the  phase  imaging. 
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Reportable  Outcomes 


We  have  published  seven  papers  in  peer-reviewed  journals  to  report  the 
research  results  from  this  project  and  presented  the  works  in  national  conferences  as 
well. 


Peer-Reviewed  Journal  Article: 

Wu  X  and  Yan  A.  2009.  Phase  retrieval  from  one  single  phase  contrast  x-ray  image, 
Optics  Express  17,  11187-11196. 

Yan  A,  Wu  X,  Liu  H.  2010.  Performance  analysis  of  the  attenuation-partition  based 
iterative  phase  retrieval  algorithm  for  in-line  phase-contrast  imaging,  Optics  Express  18: 
16074-16089. 

Yan  A  and  Wu  X.  2011 .  Apparent  linear  attenuation  coefficients  in  phase  contrast  x- 
ray  tomography,  Nuclear  Instruments  and  Methods  in  Physics  Research,  B  269  1841- 
1843. 

Yan  A,  Wu  X,  Liu  H.  201 1 .  Robustness  of  phase  retrieval  methods  in  x-ray  phase 
contrast  imaging:  A  comparison,  Medical  Physics  38:  5073-5080. 

Wong  M,  Wu  X,  Liu  H.  201 1 .  The  effects  of  x-ray  beam  hardening  on  detective 
quantum  efficiency  and  radiation  dose,  Journal  of  X-Ray  Science  and  Technology 
19:  509-519.  (In  Appendix  5) 

Wu  X,  Yan  A,  Liu  H.  2012.  X-ray  phase-shifts-based  method  of  volumetric  breast  density 
measurement,  Medical  Physics  39:  4239-4244. 

Liu  H,  Ren  Y,  Guo  H,  Xue  Y,  Xie  H,  Xiao  T,  Wu  X.  2012.  Phase  retrieval  for  hard 
X-ray  computed  tomoqraphy  of  samples  with  hybrid  compositions,  Chinese  Optics 
Letters  10:  12110  1-4. 


Published  Abstract  and  Conference  Presentation: 

X.  Wu,  A.  Yan  ,  H.  Liu:  “Contrast  mechanism  and  origin  of  artifacts  in  phase-sensitive 
x-ray  volumetric  imaging”,  Medical  Physics  36:  2785  (2009). 

X.  Wu,  A.  Yan  ,  H.  Liu:  “Contrast  mechanism  and  origin  of  artifacts  in  phase-sensitive  x- 
ray  volumetric  imaging.”  Oral  Presentation  at  the  51  st  Annual  Meeting  of  the  American 
Association  of  Physicists  in  Medicine,  July  29,  2009,  Anaheim,  California 

X.  Wu,  A.  Yan,  H.  Liu,  Improving  robustness  of  phase  retrieval  in  x-ray  phase  contrast 
imaging,  Medical  Physics  37:  3357  (2010) 

X.  Wu,  A.  Yan  ,  H.  Liu,  “Improving  robustness  of  phase  retrieval  in  x-ray  phase  contrast 
imaging,”  Oral  Presentation  at  the  52nd  Annual  Meeting  of  the  American  Association  of 
Physicists  in  Medicine,  July  19,  2010,  Philadelphia,  PA. 
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A.  Yan,  X.  Wu,  M.  A.  Yan,  X.  Wu,  and  H.  Liu:  Development  of  a  prototype  phase- 
sensitive  x-ray  breast  imaging  system”,  DoD  Breast  Cancer  Research  Program  Era  of 
Hope  2011  Conference,  Aug.  3-Aug.  5,  2011,  Orlando,  Florida. 

H.  Liu  H,  X.  Wu:  Development  and  characterization  of  phase  sensitive  x-  ray  imaging 
systems,  The  53rd  Annual  Meeting  of  the  American  Association  of  Physicists  in 
Medicine,  July  31-Aug.  4,  2012,  Vancouver,  Canada. 


Conclusion 


Aiming  to  develop  a  low-dose  and  quantitative  x-ray  phase  imaging  technique  for 
facilitating  breast  cancer  detection,  in  this  project  we  designed  and  built  a  high-voltage  x- 
ray  imaging  prototype  that  enables  x-ray  phase  imaging  and  radiation  dose  reduction,  as 
is  demonstrated  in  phantom  imaging.  We  developed  an  innovative  phase  attenuation 
duality  based  (PAD-based)  method,  which  greatly  facilitating  the  low  doe  phase  retrieval 
for  soft  tissue  body  parts  such  as  breast.  The  phantom  imaging  demonstrated  that  the 
PAD-based  phase  retrieval  greatly  enhance  the  contrast-to  noise  ratios  of  images, 
hence  it  has  good  potential  for  reducing  radiation  doses  in  breast  imaging.  We  also 
developed  a  robust  general  phase  retrieval  method,  the  attenuation-partition  based  (AP- 
based)  method  for  general  clinical  imaging.  This  method  was  demonstrated  to  be 
superior  to  the  mostly  used  phase  retrieval  method  in  literature  (TIE-based  method)  in 
performing  robust  phase  retrieval  for  low-dose  imaging.  Advancing  the  understanding  of 
the  image  contrast  mechanism  in  phase  contrast  tomography,  we  derived  a  general 
formula  relating  the  apparent  linear  attenuation  coefficients  to  the  sample's  linear 
attenuation  coefficients  and  refractive  indices  in  phase  contrast  tomography. 

Considering  the  high  breast  density  is  one  of  the  biggest  risk  factors  for  breast  cancer, 
we  proposed  for  the  first  time  an  x-ray  phase-shifts  based  method  of  breast  density 
measurement,  as  a  quantitative  imaging  method  of  the  phase-sensitive  breast  imaging. 
These  results  obtained  and  the  techniques  developed  in  this  project  contribute 
significantly  to  advancing  the  development  of  phase-sensitive  low-dose  and  quantitative 
breast  imaging. 
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Abstract:  Phase  retrieval  is  required  for  achieving  artifact-free  x-ray 

phase-sensitive  3D  imaging.  A  phase-retrieval  approach  based  on  the 
phase-attenuation  duality  with  high  energy  x-rays  can  greatly  facilitate 
for  phase  sensitive  imaging  by  allowing  phase  retrieval  from  only  one 
single  projection  image.  The  previously  derived  phase  retrieval  formula  is 
valid  only  for  small  Fresnel  propagator  phases  corresponding  to  common 
clinical  imaging  tasks.  In  this  work  we  presented  a  new  duality-based  phase 
retrieval  formula  that  can  be  applied  for  cases  with  large  Fresnel  propagator 
phases  corresponding  to  high  spatial  resolution  imaging.  The  computer 
simulation  demonstrated  superiority  of  this  new  formula  over  the  previous 
phase  retrieval  formula  in  reconstructing  the  high  frequency  components 
of  imaged  objects.  A  modified  Tikhonov  regularization  technique  has 
been  devised  for  phase  retrieval  in  cases  of  very  high  resolution  and  large 
object-detector  distance  such  that  some  Fresnel  propagator  phases  may  be 
close  or  greater  than  n.  This  new  phase  retrieval  formula  lays  the  foundation 
for  implementing  high-resolution  phase-sensitive  3D  imaging  of  soft  tissue 
objects. 

©  2009  Optical  Society  of  America 
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1.  Introduction 

Conventional  x-ray  imaging  for  soft  tissues  (such  as  breast,  brain,  liver,  etc.)  is  limited  in  its  sen¬ 
sitivity  for  detecting  subtle  tissues  pathological  changes,  since  the  imaging  relies  on  the  small 
difference  in  x-ray  attenuation  between  the  lesions  and  soft  tissues  of  variable  structure.  How¬ 
ever,  as  x-ray  wave  passes  through  tissues,  x-ray  undergoes  phase-shifts  as  well.  The  amount 
of  x-ray  phase  shift  (j)  by  tissues  is 


where  re  is  the  classical  electron  radius,  h  is  the  Plank  constant  and  c  the  speed  of  light,  E  is 
x-ray  photon  energy.  Here  pe  denotes  tissue  electron  density,  and  the  integration  of  electron 
density  pe  is  over  the  x-ray  path  and  the  integral  is  the  projected  electron  density  pe  p  [1,  2,  3]. 
It  turns  out  that  x-ray  phase-shift  differences  between  tissue  and  lesions  are  about  one  thou¬ 
sand  times  larger  than  their  attenuation  differences  [1,2,  3].  Hence  the  phase-sensitive  imaging 
has  the  potential  to  greatly  enhance  the  lesion  detection  sensitivity  and  specificity,  and  reduce 
the  radiation  dose  associated  with  x-ray  imaging.  In  the  inline  phase  contrast  imaging  as  stud¬ 
ied  in  previous  works,  the  x-rays  with  shifted  phases  diffract  from  the  object  to  the  detector, 
and  the  diffracted  x-rays  form  bright  and  dark  fringes  at  tissue  boundaries.  The  edge  enhance¬ 
ment  thus  generated  relies  on  the  spatial  coherence  of  the  x-ray  source,  and  Laplacian  and 
gradients  of  x-ray  phase-shifts  caused  by  the  object,  and  gradients  of  the  objects  attenuation 
[4,  5,  6,  7,  8,  9,  10,  11,  12,  13,  14].  The  procedure  of  disentangling  tissue  phase-shifts  from 
the  mixed  contrast  mechanism  and  retrieving  the  phase  maps  from  acquired  phase  contrast 
images  is  called  the  phase  retrieval.  Phase  retrieval  technique  plays  a  central  role  in  phase- 
sensitive  x-ray  imaging.  By  means  of  phase  retrieval,  one  can  reconstruct  a  quantitative  map 
of  phase-shifts,  a  phase  image  of  the  imaged  object  [3,  6,  13,  15,  16].  According  to  Eq.  (1), 
a  retrieved  phase  image  is  equivalently  a  map  of  imaged  objects  quantitative  projected  elec¬ 
tron  densities.  Therefore,  tissue  phase  images  can  provide  quantitative  tissue  characterization. 
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which  potentially  can  increase  the  sensitivity  and  specificity  of  a  diagnosis.  Moreover,  phase 
retrieval  is  a  necessary  procedure  for  phase-sensitive  volumetric  imaging  such  as  the  phase- 
sensitive  tomography  and  tomosynthesis  [7,  14,  15].  In  the  phase-sensitive  imaging,  x-ray  does 
NOT  propagate  along  straight  lines  because  of  x-ray  refraction  and  diffraction  at  tissue  inter¬ 
face  and  boundaries.  The  diffraction  fringes  generated  by  phase  contrast  at  tissue  interfaces 
will  obviously  render  conventional  volumetric  image  reconstruction  methods  erroneous,  since 
conventional  volumetric  reconstruction  methods  assume  that  x-rays  propagate  along  straight 
lines.  The  phase-sensitive  CT  experiments  showed  that  conventional  tomography  reconstruc¬ 
tion  algorithms  applied  to  phase-contrast  projections  can  lead  to  erroneous  tomography  images 
with  artifacts  such  as  the  streaks  and  negative  densities  [7].  The  exact  mechanism  of  the  phe¬ 
nomenon  was  revealed  in  one  of  our  recent  work.  In  fact,  if  we  denote  the  imaged  objects 
linear  attenuation  coefficient  by  p(r)  and  its  electron  density  by  pe(r),  with  the  conventional 
filtered  backprojection  method  the  reconstructed  “apparent  attenuation  coefficients”  /ireCon (/) 
was  found  as  [17]: 


(2) 


where  R2  denotes  the  object-detector  distance,  A  the  x-ray  wavelength.  Eq.  (2)  shows  that  the 
reconstructed  tomogram  consists  of  three  sub-tomograms.  The  first  sub-tomogram  is  the  tomo¬ 
gram  of  the  tissue  linear  attenuation  coefficients  p(r).  The  second  sub-tomogram  is  a  scaled 
map  of  the  3-D  Laplacian  of  tissue  electron  densities.  The  third  sub-tomogram  /jmiXed(Vj  repre¬ 
sents  the  artifacts  reflecting  the  global  variation  of  both  /r(r)  and  p,, ( r) ,  but  is  not  representing 
the  local  variations  of  ju(r)  and  pr.  ( r) .  It  is  this  term  that  causes  artifacts  such  as  the  streaks  in 
reconstruction  images.  In  order  to  correctly  deal  with  volumetric  reconstruction  issue,  as  a  gen¬ 
eral  strategy  one  should  first  disentangle  tissue  phase-shifts  from  the  mixed  contrast  mechanism 
from  acquired  phase-sensitive  projections  and  retrieve  the  phase  maps  for  all  projections,  then 
reconstruct  the  3D  images  of  object  attenuation  coefficients  and  electron  densities,  respectively. 

For  phase  retrieval  one  in  general  needs  to  acquire  at  least  two  and  mostly  four  to  five  im¬ 
ages  with  different  object-detector  distances  [3,  15].  This  requirement  is  imposed  by  the  x-ray 
wave  nature  and  its  wave  equation  [3,  8,  9,  10].  Multiple  image-acquisitions  for  a  phase  re¬ 
trieval  increase  radiation  dose  associated  with  the  technique.  For  phase  sensitive  tomography  or 
tomosynthesis,  multiple  image-acquisitions  for  each  projection-angle  make  phase-sensitive  to¬ 
mography  and  tomosynthesis  implementation  especially  cumbersome:  multiple  phase-contrast 
images  have  to  be  acquired  with  different  object-detector  distances  for  all  the  projections  dur¬ 
ing  a  scan,  and  the  number  of  required  projections  can  be  as  high  as  one  thousand.  In  addition, 
the  image  data  alignment  procedure  is  quite  time-consuming.  Moreover,  radiation  dose  is  mul¬ 
tiplied  for  multiple  image  acquisition  per  projection  view  as  well.  Therefore,  an  important  chal¬ 
lenge  is  to  find  an  effective  and  low-dose  approach  for  phase  retrievals  for  phase-sensitive  imag¬ 
ing.  Many  studies  were  conducted  addressing  this  challenge  [6,  13,  14].  We  proposed  a  unique 
solution  for  phase-sensitive  imaging  with  high  energy  x-rays.  We  noted  that  with  increasing  x- 
ray  energy  tissues  photoelectric  absorption  decreases  as  1  /En  with  rk,  3—4.  but  tissues  phase 
shift  decreases  much  slower,  and  only  as  l/E  as  is  shown  in  Eq.  (1).  Therefore,  phase-sensitive 
imaging  with  x-rays  can  fully  take  advantage  of  phase  contrast  and  greatly  reduce  radiation 
doses  associated  with  x-ray  imaging,  by  taking  advantage  of  the  dual  relation  between  the  high 
energy  x-rays  phase-shifts  and  attenuation  for  soft  tissues.  Furthermore,  analyzing  the  tissue 
attenuation  data,  we  observed  that  the  soft  tissue  attenuation  cross  sections  are  very  well  ap¬ 
proximate  by  that  of  x-ray  incoherent  scattering  for  x-rays  of  about  60  keV  to  500  keV.  Under 
this  circumstance,  we  noted  that  soft-tissue  attenuation  and  soft-tissue  phase  are  all  related  to 
the  projected  electron  density,  we  called  this  new  notion  as  the  phase-attenuation  duality.  Based 


#1 10095  -  $15.00  USD  Received  14  Apr  2009;  revised  20  May  2009;  accepted  17  Jun  2009;  published  19  Jun  2009 


Det 


S 


Fig.  1.  A  shematic  diagram  for  the  measuring  system 


on  this  concept  of  phase-attenuation  duality,  we  found  a  way  to  retrieve  a  phase-map  from  only 
one  single  recorded  image,  that  is,  the  retrieved  phase  map  is  given  by  [13]: 


(3) 


where  Ri  is  the  distance  from  the  point  x-ray  source  to  object,  AS  is  the  distance  from  object 
to  detector,  M  =  (R\  -I  Ri) / R\  is  the  geometric  magnification,  r  is  the  position  vector  (see 
Fig.  1.)  In  Eq.  3,  I(Mr\R\  +  Rt)  is  the  measured  image  intensity  at  detector  plane,  and  Im  is 
the  intensity  at  entrance,  the  operator  V2  denotes  the  Laplacian  operator.  In  Eq.  (3)  Okn  is  the 
Klein-Nishina  total  cross  section,  that  is,  the  total  cross  section  for  x-ray  Compton  scattering 
with  a  free  electron  [17].  Note  that  Okn  is  x-ray  energy  dependent,  as  is  shown  below.  We  point 
out  that  the  high  sensitivity  of  x-ray  phase  change  is  realized  as  the  high  ratio  Xrej 1 
in  this  circumstance.  In  fact,  for  clinical  imaging  photon  energy  is  no  higher  than  150  keV, 
and  in  this  case  Xrej Okn  ~  103.  In  other  words,  with  this  duality  the  high  sensitivity  of  phase 
imaging  manifests  itself  as  the  high  ratio  of  x-ray  wavelength  over  the  classical  electron  radius. 
Compared  to  the  common  multiple-acquisition  based  phase  retrieval  approaches  in  the  literature 
[3,  6,  15],  our  duality-based  phase  retrieval  method  is  a  conveniently  implemented  method 
with  tremendous  dose-saving  for  many  potential  applications  such  as  mammography.  These 
advantages  will  be  critical  for  phase-sensitive  volumetric  imaging  such  as  breast  tomosynthesis 
and  computed  tomography. 

However,  there  is  a  limitation  on  applicability  of  the  phase  retrieval  formula  Eq.  (3).  In  fact. 
Eq.  (3)  is  not  valid  for  the  cases  with  very  high  spatial  resolution  and  with  large  object-detector 
distances.  The  physics  behind  this  can  be  intuitively  explained  as  follows.  In  fact,  in  inline 
phase-sensitive  imaging  the  phase  contrast  is  formed  by  the  diffraction  from  the  object  to  the 
detector  of  x-rays  with  shifted  phases.  The  Fresnel  diffraction  can  be  thought  as  a  x-ray  wave- 
front  filtering  with  the  Fresnel  propagator  exp(inXzu2),  where  X  is  the  wavelength  and  u  is 
the  wave-front  transverse  spatial  frequency  vector,  and  z  =  Rj/M  in  terms  of  the  parameters  in 
Eq.  (3).  So  the  phase  of  the  propagator,  i.e.  nXzu2,  determines  how  strongly  the  wave-front  exit 
from  the  object  will  be  diffracted  in  x-ray  propagation.  The  phase  retrieval  formula  Eq.  (3)  was 
derived  by  linearizing  Fresnel  propagator  phase  under  the  assumption  of  nXzu 2  <C  1  [8,  13]. 
This  linearizing  is  equivalent  to  simplifying  x-ray  wave  equation  of  motion  to  the  transport  of 
intensity  of  equation  (TIE)  [8].  However,  for  cases  with  very  high  spatial  resolution  and  with 
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large  object-detector  distances,  such  as  cases  with  phase-sensitive  imaging  experiments  with 
synchrotron  based  x-ray  sources,  or  the  high  resolution  imaging  experiments  with  biological 
specimen  or  small  animals,  the  high  spatial  resolution  or  large  propagation  distances  make 
Fresnel  propagator  phases  nXzu1  larger  than  1  for  many  transverse  spatial  frequencies  of  x-ray 
wave-front.  For  example,  for  60  keV  x-rays,  if  M  =  1  and  AS  =  2m,  detector  pixel  is  5  /im,  then 
nXzu 2  =  1.299.  For  these  cases  one  cannot  linearize  Fresnel  propagators  and  should  consider 
the  full  exponential  propagators  instead.  In  section  2,  we  will  first  briefly  review  the  concept 
of  the  phase-attenuation  duality  and  how  it  greatly  simplify  the  phase  retrieval  for  soft-tissue 
objects.  We  will  then  derive  a  new  and  general  PA-Duality  based  phase  retrieval  formula,  as  the 
main  result  of  this  work.  With  this  general  formula,  one  can  retrieve  the  phase  from  one  single 
phase  contrast  image  and  is  valid  for  cases  involved  with  large  Fresnel  propagator  phases. 
Hence  this  new  phase  retrieval  formula  extends  its  applicability  to  cases  with  very  high  spatial 
resolution  and  with  large  object-detector  distances.  In  section  3,  we  show  the  simulation  results 
with  the  general  formula  for  phase  retrieval. 


2.  Phase-attenuation  duality  and  phase  retrieval  from  one  single  phase  contrast  image 
for  soft  tissue  objects 

As  is  explained  in  previous  works  [3,  4,  8,  9,  10],  in  phase-sensitive  imaging  the  imaged  object 
can  be  modeled  by  a  two-dimensional  transmission  function  T{r)  =  Ao(r)el^r\  where  0o(r) 
denotes  the  x-ray  phase-shift  caused  by  the  object,  and  Ao(?)  denotes  the  x-ray  transmission, 
or  the  attenuation-map  of  object.  For  sake  of  clarity  in  exploring  the  effects  of  non-linearity  of 
the  Fresnel  propagators,  we  assume  the  x-ray  source  is  a  quasi-monochromatic  point  source  in 
following  derivations.  We  will  discuss  more  general  x-ray  sources  in  section  3. 

Let  R\  denote  the  source-object  distance  and  AS  the  object-detector  distance.  The 
geometric  magnification  factor  M  is  equal  to  (  R\  +  AS) / R  \  .  We  also  assume  that 
one  encounters  only  moderate  variation  of  the  object  phase  function  0(r)  such  that 
|0  (r  +  XRzu/lM)  —  0  {f  —  XRju j2M)\  <C  1  [3].  Transmitting  from  the  object,  x-rays  undergo 
attenuation  and  phase-shift,  and  diffract  over  a  distance  R2  to  the  detector.  An  important  task  is 
to  find  out  the  formula  of  the  image  intensity  as  function  of  the  objects  transmission  function 
and  imaging  geometry.  A  general  solution  had  been  found  in  one  of  our  previous  works  [3]. 
Starting  from  the  paraxial  Fresnel-Kirchhoff  diffraction  theory,  we  found  out  in  this  work  that 
the  Fourier  transform  of  x-ray  irradiance  at  detector  plane  is  given  by: 
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where  I(jd)  denotes  the  x-ray  irradiance  at  the  detector  entrance,  and  the  symbol  -IX  (■)  denotes 
the  2-D  Fourier  transform  and  Im  is  the  intensity  of  the  incident  x-ray  upon  the  object,  and 
u  is  the  spatial  frequency  vector  in  the  object  plane.  Note  that  the  Fresnel  propagator  phases 
nXRiu2 /M  enter  Eq.  (4)  as  the  argument  of  sinusoidal  transfer  functions.  For  cases  with  mod¬ 
erate  spatial  resolutions  u  and  AS  ,  Eq.  (4)  is  much  simplified,  and  the  inverse  Fourier  transfor¬ 
mation  of  the  simplified  Eq.  (4)  can  be  written  as  [8]: 


I(r-Rl+R2) 


2  kM 


(5) 


Equation  (5)  is  the  TIE-based  formula  for  image  irradiance,  which  was  first  derived  by  Nugent 
and  colleagues [3].  Recently  Guigay  et  al  showed  the  limitations  of  Eq.  (5)  with  theoretical 
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analysis  and  experiments  [17].  They  found  that  for  in-line  phase-sensitive  imaging  with  long 
R2  and  high  resolutions  such  that  nXR2u2  /M  >  1  ,  the  phase  retrieval  based  on  Eq.  (5)  would  be 
inaccurate,  one  has  to  use  the  general  formula  Eq.  (4)  for  accurate  quantitative  phase  imaging. 

Furthermore,  we  note  that  soft  tissues  encountered  in  clinical  imaging  are  composed  almost 
exclusively  by  atoms  of  the  light  elements  with  the  atomic  numbers  Z  <  10.  For  example,  the 
composition  of  breast  tissue  consists  mainly  of  light  elements  such  as  hydrogen,  carbon,  nitro¬ 
gen  and  oxygen.  In  contrast,  the  sum  of  weight  fractions  of  other  heavier  elements  in  breast  tis¬ 
sue  is  only  less  than  1%  [18].  The  same  is  true  for  other  soft  tissues  such  as  the  brain  grey/white 
matter  [18].  Analyzing  the  tissue  attenuation  data,  we  observed  that  the  soft  tissue  attenuation 
cross  sections  are  very  well  approximated  by  that  of  x-ray  incoherent  scattering  for  x-rays  of 
about  60keV  to  500keV[19].  In  addition,  the  incoherent  scattering  functions  become  linearized 
as  well  [19].  Under  this  circumstance,  while  the  phase-shift  is  caused  by  the  coherent  x-ray 
scattering,  but  soft  tissue  attenuation  is  determined  by  the  incoherent  x-ray  scattering  for  soft 
tissues,  so  both  the  x-ray  phase-shift  and  attenuation  by  soft  tissues  are  all  determined  by  the 
projected  electron  density  for  these  high  energy  x-rays.  We  called  this  complementary  relation¬ 
ship  between  phase-shift  and  attenuation  for  soft  tissues  as  the  phase-attenuation  duality.  When 
the  P-A  duality  holds,  x-ray  attenuation  and  phase  shift  by  the  object  is  related  in  following 
way: 

Hr)  =  -ArfpejP(r),  A20(r)  =  exp(-OkNpe,„(r)),  (6) 

where  Okn  is  the  total  cross  section  for  x-ray  Compton  scattering  with  a  free  electron: 
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where  rj  =  Ephoton/mec2,  here  we  denote  the  photon  energy  of  the  primary  x-ray  beam  by 
/-photon,  and  mec2  is  the  rest  electron  energy  and  equal  to  51  IkeV.  The  duality  condition  can  be 
used  to  simplify  Eq.  (4).  In  fact,  we  can  rewrite  Eq.  (4)  as 
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Especially  with  the  duality  condition  we  found 
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In  Eq.  (9)  we  made  an  approximation  of  dropping  the  negligible  contribution  of  the  higher 
order  terms  of  Okn/^tv  The  approximation  is  based  on  GknMt>  ~  10  and  we  assumed 
|p  (r  +  XR2u/2M)  —  0  (r  —  XR2u/2M)\  -C  1,  moderate  variations  of  object  phase  shifts.  Sub¬ 
stituting  Eqs.  (9)  into  Eq.  (8),  we  found: 
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Equation  (10)  is  the  central  result  of  this  work,  it  is  the  formula  for  x-ray  irradiance  for  high 
resolution  soft  tissue  imaging  with  high  energy  x-rays.  In  the  derivation  the  parameter  A re! Okn 
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plays  a  significant  role.  It  is  interesting  to  note  that  Gureyev  et  al  obtained  a  similar  result  in 
ref.[20]  (Eq.(24))  with  the  exception  that  the  term  nXRiu1  /  M  in  Eq.  (10)  was  absent.  The  dif¬ 
ferences  may  lie  in  the  different  assumptions  made  regarding  the  object  transmission  functions 
and  x-ray  energies.  Gureyev  and  coauthors  assumed  the  objects  as  single  material  homoge¬ 
neous  objects  [20],  but  our  results  are  for  heterogeneous  objects  of  the  light  elements  imaged 
with  high  energy  x-rays  (60  keV  or  higher),  i.e.,  for  the  cases  the  phase-attenuation  duality 
holds.  In  addition,  in  their  derivation  they  separated  the  slowly  and  rapidly  varying  components 
of  the  amplitude  and  phase,  and  they  assumed  that  the  magnitude  of  the  rapidly  varying  com¬ 
ponents  are  much  smaller  than  1,  but  we  did  not  make  such  limiting  assumption  in  derivation  of 
Eq.  (10).  We  stress  that  the  Klein-Nishina  cross-section  is  energy  dependent  shown  in  Eq.  (7) 
and  X  re  j (Jkn  is  about  1000  for  inline  x-ray  energy  range  used  in  diagnostic  imaging.  This  large 
Are/<TKN  represents  the  high  sensitivity  of  phase  contrast  in  tissue  differentiation,  but  it  also 
makes  the  x-ray  irradiance  formula  Eq.  (10)  getting  very  much  simplified.  It  is  easy  to  verify 
that  for  the  low  resolution  cases  such  as  those  with  n'kR^u1  jM  <C  1,  Eq.  (10)  is  reduced  to 


we  recover  our  previous  result  as  shown  in  [13].  On  the  other  hand,  for  general  cases 
nXR^u2  /M  ~  1  or  nXRiu1  /M  >  1  Eq.  (11)  cannot  be  applied,  and  one  should  use  Eq.  (10) 
for  phase  retrieval  as  is  shown  below. 

3.  Simulation  results 

Comparing  Eq.  (10)  and  Eq.  (11)  for  use  in  performing  phase  retrieval,  the  only  difference  lies 
in  the  inverse  filters  derived  for  retrieving  JF  (Aq).  For  phase  retrieval  based  on  Eq.  (10)  the 
inverse  filter  is  found  as 


1  _  1 
cos(a)  +  (2y+a)sin(a)  <5-cos(ai)’ 


(12) 


where  a  =  tiXRju1  j  M ,  the  FT-space  Fresnel  Propagator  phase  as  defined  in  section  1,  and 
y  =  X  re  / <7k x •  which  represents  the  phase  contrast  enhancement  as  explained  earlier.  8  = 
1  T  (2y+  a)2  and  ai  =  a  —  sin-1  [(2y+  a) /8\.  On  the  other  hand,  for  phase  retrieval  based 
on  Eq.  (11)  the  inverse  filter  will  be 
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Since  y  ~  103,  the  filters  in  Eq.  (12)  and  Eq.  (13)  are  close  when  a  is  small.  When  a  is  getting 
bigger,  close  to  or  greater  than  n  (oq  close  to  n/2),  the  difference  between  the  two  inverse 
filters  is  getting  bigger.  Figure  2  shows  a  plot  of  the  two  filters,  with  respect  to  the  change  of  a. 
The  inverse  filter  in  Eq.  (12),  depicted  in  Fig.  2  as  the  solid  line,  changes  the  sign  when  a  passes 
through  kn  (and  thus  a.\  passes  through  (k—  \/2)n)  and  diverges  to  °°  when  ai  =  (k—  1/2 )n, 
k=  1,2,  •  •  • .  On  the  other  hand,  the  filter  in  Eq.  (13),  depicted  as  the  dashed  line  in  Fig.  2,  keeps 
positive  and  approaches  zero  as  a  gets  bigger.  So  the  difference  of  the  retrieved  phase  maps 
using  Eq.  (10)  and  Eq.  (11)  will  be  highly  dependent  on  the  value  of  the  Fresnel  propagator 
phase  a  =  nXR2u2/M. 

To  demonstrate  this  phenomenon,  we  performed  two  simulations.  First,  we  set  the  phase  shift 
(j>  to  be  a  sine-function  in  x  direction  with  period  of  32  pixel  size  and  constant  in  y  direction. 
Assuming  the  pixel  pitch  is  p  jum,  the  spatial  frequency  corresponding  to  0,  in  x  direction,  is 
RrO  =  l/32y>  1/jum.  In  our  simulations,  for  the  purpose  of  simplification,  the  magnification  M 
as  well  as  the  intensity  at  entrance,  /m,  are  assumed  to  be  1.  The  x-ray  photon  energy  is  set 
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Fig.  2.  Plots  of  filters  (13)  and  (12) 


to  60  keV.  This  gives  the  x-ray  wavelength  A  =  2.0667  x  10  11  m  and  the  total  cross  section 
for  x-ray  Compton  scattering  with  a  free  electron  <7kn  =  5.4526  x  1 0  29  m2  as  determined 
from  Eq.  (7).  The  classical  electron  radius  re  is  set  to  2.82  x  10~15  m.  By  assuming  the  de¬ 
tectors  resolution  p  =  1  /im  and  the  phase  shift  <l>(x.y)  =  0.5sin(32/?x)  —  1,  according  to  the 
phase-attenuation  duality  discussed  earlier,  the  attenuation  map  is  related  to  the  phase  map 
as  A/(x,y)  =  exp  [OYm/Xre  •</>]  =  exp  [(sin(32px)  —  2) /2y],  with  y=  Aty/okN  =  1-0726  x  103. 
We  then  performed  Fresnel  diffraction  to  determine  the  phase-contrast  image  irradiances  at  dis¬ 
tances  downstream [21],  and  the  retrieved  maps  of  4/  were  found  by  using  the  inverse  filters  of 
Eq.  (12)  and  Eq.  (13),  which  were  derived  based  on  Eq.  (10)  and  Eq.  (11),  respectively.  Because 
of  the  large  value  of  y,  phase  contrast  image  irradiances  {/}  were  found  be  maps  with  each  I  be¬ 
ing  dominated  by  two  single  frequency  component  at  ux  =  ±;/aq  (uxq  =  1  /32 p  =  1/32  x  106  m) 
plus  a  zero  frequency  component.  This  being  so,  the  shape  and  variation  of  the  retrieved  Ag- 
maps  were  determined  by  their  dominant  frequency  component  at  frequencies  ux  =  ±ux 0.  In  or¬ 
der  to  clearly  show  the  differences  between  phase  retrievals  based  on  Eq.  (10)  and  Eq.  (1 1),  we 
examined  the  ratio  of  the  dominant  frequency  components  of  the  two  retrieved  Ajj-maps  based 
on  Eq.  (10)  and  Eq.  (11),  respectively.  We  found  the  ratio  being  r  =  (1  +2ya)/(<5  -cos(ai)) 
for  a  given  object-detector  distance  Ri-  This  ratio  reflects  the  capability  of  achieving  accu¬ 
rate  phase  retrieval  with  Eq.  (10)  against  Eq.  (11).  For  small  1,  that  is,  for  small  Fres¬ 
nel  propagator  phases,  this  ratio  is  1,  so  Eq.  (10)  against  Eq.  (11)  both  are  equally  good  for 
the  accurate  phase  retrieval.  However,  when  a  is  chosen  such  that  cosjcq)  is  close  to  0,  this 
ratio  can  be  very  large.  For  example  by  setting  R2  =  49.55  m,  we  have  a  mentioned  above 
equals  to  n,  the  ratio  r  =  (1  +2y)/(t>  ■  cos(ai))  ss  6700.  This  large  ratio  means  that  the  re¬ 
trieved  ,4 /-map  based  on  Eq.  (11)  has  diminished  image-contrast  and  erronous  values,  as  is 
shown  in  Fig. 3(c).  The  phase  retrieval  result  based  on  Eq.  (10),  shown  in  Fig.  3(d),  how¬ 
ever  is  much  more  accurate.  When  R2  =  49.55  m,  the  FT-space  Fresnel  Propagator  phase 
a  =  nXRyu2 /M  ss  n.  So  we  have,  from  Eq.  (10),  $  ((M2//in)/)  ss  —  &  (Aq),  which  is  differ¬ 
ent  from  the  result  &  ((M2 //;„)/)  «  [1  +27r(Are/oKN)]  &  (Ag)  derived  from  the  conventional 
TIE  prediction  (Eq.  (11)),  which  does  not  hold  for  such  large  Fresnel  Propagator  phases.  This 
phenomenon  is  shown  in  Fig.  (3)(b),  where  the  intensity  I  is  acquired  with  Fresnel  diffraction, 
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in  which  the  dark-white  area  in  A5  is  reversed.  Note  that  the  inverse  filter  in  Eq.  (12)  is  singular 
at  oq  =  [k—  1/2)  tc,  k  =  1 ,2,  •  •  • ,  because  of  the  zero-crossing  of  cos(ai).  To  deal  with  the  sin¬ 
gularity  associated  with  the  zero-crossing  of  cos(ai),  one  should  regularize  the  inverse  filter. 
A  common  regularization  technique  is  the  Tikhonov  regularization,  which  modifies  1  / cos(ai) 
in  the  filter  of  Eq.  (12)  to  cos(ai)/(cos(ai)2  +  jc2),  where  K  is  the  Tikhonov  regularization 
parameter.  In  this  specific  example  of  the  sinusoidal  object,  however,  since  we  know  &  (/)  is 
dominated  by  frequency-components  only  at  ux  =  ±ux o  and  ( u ,  v)  =  (0, 0),  by  setting  &  (/)  =  0 
at  all  other  frequency  points  purposely,  we  can  avoid  the  problem  caused  by  the  zero-crossing 
of  the  inverse  filter. 


(a)  (b)  (c)  (d) 

Fig.  3.  Comparation  of  the  quality  of  the  retrieved  results  Aq  using  Eq.  (11)  and  Eq.  (10). 
A5  =  exp[0/y],  where  phase  shift  (j)  =  0.5sin(32px)  —  1,  pixel  pitch  p=  1pm,  and  object  to 
detector  distance  R2  =  49.55m.  (a)  is  the  true  Ajp  (b)  is  the  phase-contrast  intensity  image 
I  acquired  from  Fresnel  diffraction;  (c)  is  retrieved  A(2  using  Eq.  (11);  (d)  is  retrieved 
using  Eq.  (10). 


In  another  simulation,  we  set  all  the  parameters  the  same  as  the  previous  simulation  except 
replacing  the  test-object  with  the  well-known  Lena  image  as  is  shown  in  Fig.  4(a).  The  simula¬ 
tion  results  for  the  case  of  R2  =  5  m  are  shown  in  Figures  4(c)-4(d).  The  phase-contrast  image 
I,  shown  in  Fig.  4(b),  is  acquired  with  Fresnel  diffraction  from  the  object  downstream  to  the 
detector  [21].  Figure  4(c)  shows  the  retrieved  image  with  the  inverse  filter  of  Eq.  (13),  which 
is  based  on  Eq.  (11).  Since  the  filter  of  Eq.  (13)  is  incorrect  when  the  Fresnel  propagator  phase 
a  is  big,  say  close  or  greater  than  n,  the  high  frequency  components  of  the  object  was  unduly 
suppressed  in  the  phase-retrieval.  This  deficiency  is  reflected  in  Fig.  4(c)  such  that  the  retrieved 
image  gets  blurred  in  locations  with  high  contrast  changes,  such  as  the  hair,  the  edge  of  the  hat 
and  body.  In  contrast,  the  phase  retrieval  using  the  inverse  filter  of  Eq.  (12),  which  is  based  on 
Eq.  (10),  preserves  objects  high  frequency  components,  as  is  shown  in  Fig.  4(d). 

As  stated  above,  a  regularization  must  be  used  to  overcome  the  difficulty  associated  with  the 
zero-crossing  in  the  denominator  of  inverse  filter  in  Eq.  (13).  The  zero-crossing  occurs  obvi¬ 
ously  at  a\  =  (k  1  /2) 7T,  k  =  1 , 2,  •  •  • .  The  Tikhonov  regularization  would  replace  1  / cos(ai ) 
in  the  inverse  filter  of  Eq.  (13)  by  the  following 


cos(ai) 
cos(a2)2  +  k 2  ’ 


(14) 


where  K  is  a  small  constant.  However,  this  form  of  regularization  is  not  optimal  since  it  will 
deviate  from  the  original  inverse  filter  when  oci  is  small.  In  other  words,  the  test-object’s  low 
frequency  components  would  not  be  fully  reconstructed  in  the  retrieval.  In  order  to  get  around 
this  drawback,  we  note  that  the  zero-crossing  occurs  only  when  oq  is  close  or  greater  than 
7r/2,  while  the  most  low  frequency  components  have  a  <  1.  Hence  we  set  the  regularization 
parameter  K  =  0  when  a  is  less  than  some  threshold  value  Oq,  and  K  >  0  when  a  is  greater  than 
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Fig.  4.  Comparation  of  the  retrieved  results  using  Eq.  (11)  and  Eq.  (10)  with  Ri  =  5m.  (a) 
is  the  true  Aq;  (b)  is  the  phase-contrast  image  1  acquired  from  Fresnel  diffraction;  (c)  is 
retrieved  using  Eq.  (11);  (d)  is  retrieved  using  Eq.  (10). 


Oq.  In  our  simulation,  we  set  this  threshold  ocq  =  0.8,  so  K  =  0  when  a  <  0.8,  and  jc  =  0.12  when 
a  >  0.8.  We  can  see  from  the  image  in  Fig.  4(c)  that  most  of  the  low  frequency  information 
is  preserved.  On  the  other  hand,  the  quality  of  the  results  can  also  be  numerically  determined 
by  computing  the  error  comparing  to  the  true  Aq  image.  By  taking  the  /^-measure,  which  is 
defined  as  the  integrated  squared  errors,  we  have  the  absolute  errors  for  methods  (10)  and  (11) 
being  0.014  and  0.021  respectively. 

4.  Conclusions 

Phase  retrieval  is  required  for  achieving  artifact-free  x-ray  phase-sensitive  3D  imaging.  Multi¬ 
ple  image-acquisitions  per  projection  are  required  in  general  for  phase  retrieval  [3,  4,  15,  17, 
22,  23,  24].  This  requirement  impedes  the  implement  of  x-ray  phase-sensitive  3D  imaging, 
especially  for  phase-sensitive  3D  imaging  in  clinical  applications.  A  phase-retrieval  approach 
based  on  the  phase-attenuation  duality  with  high  energy  x-rays  can  greatly  facilitate  for  phase 
sensitive  imaging  by  allowing  phase  retrieval  from  only  one  single  projection  image.  However, 
the  previously  derived  phase  retrieval  formula  is  valid  only  for  small  Fresnel  propagator  phases 
corresponding  to  common  clinical  imaging  tasks.  In  this  work  we  generalized  the  previous 
phase  retrieval  formula  to  Eq.  (10),  which  can  be  applied  for  cases  with  large  Fresnel  propaga¬ 
tor  phases  corresponding  to  high  spatial  resolution  imaging.  The  computer  simulation  demon¬ 
strated  superiority  of  this  new  formula  over  the  previous  duality-based  phase-retrieval  formula 
in  reconstructing  the  high  frequency  components  of  imaged  objects.  In  addition,  a  modified 
Tikhonov  regularization  technique  has  been  devised  for  phase  retrieval  in  cases  of  very  high 
resolution  and  large  object-detector  distance  such  that  some  Fresnel  propagator  phases  may  be 
close  or  greater  than  n.  This  new  phase  retrieval  formula  of  Eq.  (10)  lays  the  foundation  for 
implementing  high-resolution  phase-sensitive  3D  imaging  of  soft  tissue  objects. 
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Abstract:  The  phase  retrieval  is  an  important  task  in  x-ray  phase  contrast 
imaging.  The  robustness  of  phase  retrieval  is  especially  important  for 
potential  medical  imaging  applications  such  as  phase  contrast  mammogra¬ 
phy.  Recently  the  authors  developed  an  iterative  phase  retrieval  algorithm, 
the  attenuation-partition  based  algorithm,  for  the  phase  retrieval  in  inline 
phase-contrast  imaging  [1].  Applied  to  experimental  images,  the  algorithm 
was  proven  to  be  fast  and  robust.  However,  a  quantitative  analysis  of  the  per¬ 
formance  of  this  new  algorithm  is  desirable.  In  this  work,  we  systematically 
compared  the  performance  of  this  algorithm  with  other  two  widely  used 
phase  retrieval  algorithms,  namely  the  Gerchberg-Saxton  (GS)  algorithm 
and  the  Transport  of  Intensity  Equation  (TIE)  algorithm.  The  systematical 
comparison  is  conducted  by  analyzing  phase  retrieval  performances  with 
a  digital  breast  specimen  model.  We  show  that  the  proposed  algorithm 
converges  faster  than  the  GS  algorithm  in  the  Fresnel  diffraction  regime, 
and  is  more  robust  against  image  noise  than  the  TIE  algorithm.  These 
results  suggest  the  significance  of  the  proposed  algorithm  for  future  medical 
applications  with  the  x-ray  phase  contrast  imaging  technique. 
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1.  Introduction 

Differing  from  the  conventional  x-ray  imaging,  which  relies  on  the  small  differences  in  x-ray 
attenuation  changes  between  tissues  variable  structure,  inline  phase  contrast  imaging  is  based 
on  the  tissues’  phase-shifts  diffraction  from  the  object  to  the  detector.  Since  x-ray  phase-shift 
differences  between  tissue  and  lesions  are  about  one  thousand  times  larger  than  attenuation 
differences  [2,  3,  4],  x-ray  phase  contrast  imaging  has  the  potential  to  enhance  the  lesion  detec¬ 
tion  sensitivity  and  specificity,  and  reduce  the  radiation  dose  compared  to  conventional  x-ray 
imaging.  In  the  inline  phase  contrast  imaging,  the  diffracted  phase-shifts  form  bright  and  dark 
fringes  at  tissue  boundaries  and  this  bright  and  dark  fringes  are  called  edge  enhancement.  The 
edge  enhancement  relies  on  the  spatial  coherence  of  the  x-ray  source,  the  Laplacian  and  gra¬ 
dients  of  x-ray  phase-shifts  caused  by  the  tissue,  and  the  gradients  of  the  objects  attenuation 
[5,  6,  7,  8,  9,  10,  11,  12,  13,  14,  15],  One  procedure  of  phase  contrast  imaging  is  to  disentangle 
tissue  phase-shifts  from  the  mixed  contrast  mechanism  and  recover  the  phase  maps  from  ac¬ 
quired  phase  contrast  images.  This  procedure  is  called  phase  retrieval.  Phase  retrieval  technique 
plays  a  central  role  in  phase  contrast  x-ray  imaging.  By  means  of  phase  retrieval,  one  can  recon¬ 
struct  a  quantitative  map  of  phase-shifts,  a  phase  image  of  the  imaged  object  [4,  7,  14,  16,  17]. 
The  amount  of  x-ray  phase-shifts  (j)  by  tissues  is  determined  by 

0(?)  =  -  (jf)  re  J Pe(r,z)dz=-  rePe,p(r),  (1) 

where  re  is  the  classical  electron  radius,  h  the  Plank  constant,  c  the  speed  of  light,  E  the  x-ray 
photon  energy,  and  pe  p,  the  integration  of  the  electron  density  pe  over  the  x-ray  path,  is  called 
the  projected  electron  density  [2,  3, 4].  So  a  retrieved  phase  map  is  equivalently  a  map  of  imaged 
object’s  quantitative  projected  electron  densities.  Moreover,  phase  retrieval  is  also  a  necessary 
procedure  for  phase-sensitive  volumetric  imaging,  such  as  the  phase  sensitive  tomography  and 
tomosynthesis,  to  acquire  the  artifact  free  3D  images  of  object  attenuation  coefficients  and 
electron  densities  [8,  15,  16]. 

Phase  retrieval  is  based  on  x-ray  propagation  equation  derived  either  from  the  Fresnel  diffrac¬ 
tion  or  the  Wigner  distribution  based  phase-space  formalism  [5,  18,  9,  19,  20].  To  be  specific, 
let  0(r)  be  the  x-ray  phase-shift  caused  by  the  imaged  object,  and  Ao(r)  the  x-ray  transmis¬ 
sion,  or  the  attenuation-map  of  the  object.  Then  the  Fourier  transform  of  the  x-ray  intensity, 
JF  (/)  ( u )  =  fR 2  7(3 c)  exp  [27 lix-  u]  dx,  at  position  away  from  the  object  with  distance  7?2,  of  the 
monochromatic  point  x-ray  source  starting  at  a  place  away  from  the  object  R\  distance,  can  be 
modeled  by  the  following  [9] 
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where  /m  is  the  incident  x-ray  intensity  at  R\.  X  the  wavelength  of  the  monochromatic  point  x- 
ray  source  and  M  =  (R\  +  Rj) / R\  is  the  geometric  magnification.  When  the  FT-Space  Fresnel 
propagator  kXRju2  /M  <C  1,  Eq.  (2)  can  be  simplified  to  the  Transport  of  Intensity  Equation 
(TIE)  [21,4,  9] 

/(M,  +*)  =  ^  {aJ  (£)  -  ^V.  (l)  J .  (3) 

It  is  worthy  to  note  that  Eq.  (3)  is  valid  only  for  moderate  resolution  images.  For  high  resolution 
images,  i.e.  when  the  FT-Space  Fresnel  propagator  nXRiu2 /M  is  close  or  greater  than  n,  any 
phase  retrieval  algorithms  based  on  Eq.  (3)  need  to  be  substituted  to  Eq.  (2)  [22,  23,  24].  In  this 
paper,  the  algorithms  discussed  are  all  based  on  Eq.  (3),  i.e.  the  moderate  resolution  is  satisfied. 

Numerous  algorithms  have  been  suggested  on  how  to  effectively  recover  the  phase-shift  from 
the  phase  contrast  images.  Among  these,  two  algorithms  are  most  widely  used.  One  is  the  TIE 
algorithm  implemented  by  Allen  and  Oxley  in  [25].  The  other  is  the  GS  algorithm  developed 
first  by  Gerchberg  and  Saxton  in  [26]  and  later  improved  by  Fienup  [27,  28].  These  two  algo¬ 
rithms  have  both  their  advantages  and  disadvantages.  The  TIE  algorithm  is  a  direct  approximate 
method  which  is  fast  but  is  unstable  with  noisy  data;  the  GS  algorithm  on  the  other  hand  is  rel¬ 
atively  more  stable  than  the  TIE  algorithm[25]  but  the  converging  speed  is  slow  especially  for 
the  field  of  medical  imaging.  In  [1],  the  authors  developed  an  Attenuation-Partition  Based  Al¬ 
gorithm  (APBA)  based  on  the  phase-attenuation  duality  property  of  soft  tissues  under  higher 
x-ray  energy.  This  algorithm  is  fast  and  stable  for  potential  medical  imaging.  We  compared  the 
performance  of  this  algorithm  with  the  TIE  algorithm  for  two  groups  of  data  under  the  condi¬ 
tion  of  medical  imaging  in  [1],  including  the  phase  retrieval  from  phase-contrast  images  of  a 
breast  lumpectomy  specimen.  In  this  paper,  we  will  make  a  systematic  analysis  about  this  al¬ 
gorithm  and  compare  its  performance  with  the  one  of  the  GS  algorithm  and  the  TIE  algorithm 
with  simulated  data. 

The  paper  is  organized  as  follows.  In  Section  2,  we  first  summarize  the  attenuation-partition 
based  algorithm  (APBA),  which  is  motivated  by  our  observation  of  the  phase-attenuation 
duality [14].  Then  we  give  a  measure,  called  total  variation,  used  to  evaluate  the  closeness  of 
two  image  data.  This  measure  is  used  as  a  quantitative  standard  in  comparing  the  performance 
between  different  algorithms  in  the  following  section.  In  Section  3,  we  first  develop  a  breast 
specimen  model  which  can  reflect  the  attenuation  and  phase  changes  with  respect  to  the  x-ray 
energy  change  (Section3.1),  and  then  compare  the  performance  of  the  algorithm  with  the  GS 
algorithm  (Section  3.2)  and  the  TIE  algorithm  (Section  3.3).  Finally,  we  conclude  this  paper  in 
Section  4. 


2.  The  attenuation-partition  based  algorithm  (APBA)  and  an  image  accuracy  measure 

2.1.  The  Attenuation-Partition  Based  Algorithm 

The  attenuation-partition  based  algorithm  (APBA)  is  a  recently  developed  iterative  algo¬ 
rithm  for  phase  retrieval [  1],  It  was  derived  from  our  previous  notion  of  the  phase-attenuation 
duality[14],  and  it  takes  advantage  of  the  correlation  between  the  attenuation  and  phase-shift  for 
phase  retrieval.  As  is  well  known,  tissue’s  attenuation  change  Ao  in  the  diagnostic  x-ray  imag¬ 
ing  arose  from  three  x-ray  and  tissue  interactions:  the  photoelectric  absorption,  the  coherent 
scattering,  and  the  incoherent  scattering [29,  30,  9,  14].  However,  among  the  three  interactions, 
the  attenuation  caused  by  incoherent  scattering  Akn,  which  is  dominated  by  the  x-ray  Compton 
scattering,  deserves  a  special  attention.  This  is  because  both  of  Akn  and  the  x-ray  phase-shift  <j> 
are  determined  by  the  tissue’s  projected  electron  density; 


Akn(?)  =  exp 


0(c)  —  Arepep, 


(4) 
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where  A  is  the  x-ray  wavelength,  re  the  classical  electron  radius,  pe  p  the  projected  electron 
density  as  defined  in  Eq.  (1),  and  o^n  is  the  total  cross  section  for  x-ray  Compton  scattering 
with  a  free  electron: 


t^KN  (Ephoton  )  —  A?Tre 


1+1  [2(i+i) 

tj2  [  l+2rj  77 


+  —  log(l  +277) 


1  +  377  1 
(I+277)2/’ 


(5) 


with  rj  =  Ephoton/mec2.  Here  we  denote  the  photon  energy  of  the  primary  x-ray  beam  by  Ephoton 
and  mec2  is  the  rest  electron  energy.  Eq.(4)  suggests  clearly  that  the  x-ray  attenuation  and  phase- 
shift  by  tissue  has  certain  correlation.  Our  idea  is  to  utilize  this  correlation  for  facilitating  the 
phase  retrieval. 

Of  course,  the  extent  of  this  correlation  between  phase  and  attenuation  depends  on  the  x- 
rays  photon  energy  as  well  as  the  tissues  physical  composition.  For  example,  for  light  elements 
with  atomic  numbers  Z  <  10,  x-ray  attenuation  is  dominated  by  the  Compton  scattering  for 
x-rays  of  60  keV  or  higher,  i.e.  +0  ~  ,4  kn  [  1 4] .  We  call  this  relationship  between  phase-shift 
and  attenuation  the  phase-attenuation  duality.  The  phase-attenuation  duality  can  be  used  for 
phase  retrieval  as  follows.  Consider  a  phase  contrast  imaging  setting  with  a  point  source  of 
wavelength  A.  The  object  is  at  a  distance  R\  from  the  source.  We  denote  AS  as  the  distance  from 
object  to  detector  plane,  M  =  (R\  +  R2) /R 1  the  geometric  magnification  factor,  Im  the  entrance 
x-ray  intensity  at  R\ ,  and  Id(jd)  the  x-ray  intensity  at  the  detector  plane.  For  convenience, 
we  denote  I  =  M 2  ■  Iq (jD)/kn  as  the  normalized  intensity  of  phase-contrast  image.  When  the 
phase-attenuation  duality  holds,  the  phase  map  (j)  (r)  can  be  robustly  retrieved  from  just  a  single 
projection  image  [14]: 


A|N(r)  =  S)(7)  =  #-1 


(  ^00  \ 

y  1  +  47T2Au2  J  ’ 


(6) 


where 

l=^Rb-—'  (7) 

2nM  Okn 

and  A),  for  sake  of  convenience,  is  called  the  “duality  transform”  acting  on  the  normalized 
image  I. 

In  general  imaging  cases,  such  as  with  low  energy  x-rays  or  an  object  contains  calcified 
tissues  such  as  calcification,  this  phase-attenuation  duality  does  not  hold.  However,  we  can  still 
factor  out  tissue’s  total  attenuation  Ao  as 


A0(?)  =  Akn(?)  •Ape)Coh(r),  (8) 

where  we  denote  the  attenuation  caused  by  photoelectric  absorption  and  coherent  scattering 
by  APe,coh-  Strictly  speaking,  Okn  is  only  Compton  scattering  cross-section,  it  may  be  slightly 
different  from  the  incoherent  scattering  cross-section  for  high-Z  elements.  This  is  because  while 
Compton  scattering  is  x-ray  scattering  from  the  free  electrons,  the  incoherent  scattering  is  that 
from  the  bound  atomic  electrons [31].  So  when  we  factor  Ao  =  Akn  'Ape.coh  (Eq.  (8))  we  actually 
factor  the  small  residual  binding  effect  of  atomic  electrons  into  ApeiCOh-  With  this  understanding, 
Eqs.(4)  and  (8)  are  rigorously  valid.  The  notion  of  Eq.  (6)  and  Eq.  (8)  led  us  to  the  development 
of  the  attenuation-partition  based  algorithm  [1],  While  the  derivations  and  the  algorithm  details 
of  this  method  can  be  found  from  Ref.  [1],  a  brief  description  of  the  method  is  as  follows. 
Our  goal  is  to  retrieve  the  phase  map  from  the  two  normalized  images:  one  is  the  object’s 
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attenuation  image  A5  acquired  with  Ri  =  0,  and  the  other  is  the  acquired  phase  contrast  image 
I  =  M 2  ■  //j(r'/j)  //in  with  Ri  >  0.  Employing  the  acquired  phase  contrast  image  /  and  the  duality 
transform  Eq.  (6),  we  will  first  obtain  an  estimate  for  the  attenuation-component  Akn  and  phase 
map  (j).  We  then  rewrite  Eq.  (8)  as 

^0=^KN  —  8a,  8A  =  Akn(1  —  Ape.coh),  (9) 

and  find  the  correction  term  (5,4  using  the  estimate  of  Akn-  We  then  employing  the  Fresnel 
Diffraction  transform  (as  defined  in  Eq.  (11))  to  transport  the  wavefield  8Ae'^  from  /\’|  to  AS. 
We  can  find  81  =  (8Ael&)  |2,  which  is  the  difference  between  phase  contrast  image  /  and 

the  “duality-only”  counterpart  /kn  =  |3r  (Akn<?!<^)|~-  With  the  corrected  “duality-only”  im¬ 
age  intensity  /kn  =  [\Tl  +  \/8 /  j  we  can  start  a  new  round  of  iterations  by  repeating  above 
procedure.  For  a  rigorous  analysis  of  the  iterative  algorithm  and  its  convergence  interesting 
readers  are  referred  to  [1].  Note  that  the  equation  \/l  =  \JIkn  —  \/<5/  is  generally  valid,  since 
it  is  actually  a  result  of  x-ray  Fresnel  diffraction  and  extremely  smallness  of  hard  x-ray  wave¬ 
length  compared  to  finest  resolution  achievable  in  the  phase  contrast  imaging.  While  interesting 
readers  can  find  the  mathematical  proof  of  this  equation  in  [1],  an  intuitive  explanation  of  this 
formula  comes  from  the  x-ray  propagation.  In  such  a  wave  propagation  process,  or  the  so-called 
semiclassical  wave  propagation,  the  phase  of  a  wave  field  evolves  essentially  according  to  the 
free-space  Hamilton- Jacobi  equation  from  its  initial  phase  value.  So  if  we  denote  the  Fresnel 
free  propagation  as  a  Fresnel  transform  gr  acting  on  the  initial  wavefield,  therefore  the  two 
resulted  wavefields  (Akn  exp [;</)] )  and  (5Aexp[i0])  would  have  the  same  resultant  phases, 
so  the  resultant  amplitude  from  the  two-wave  superposition  is  given  as  \Jl  =  \/Ikn  —  Vo7. 

The  above  iterative  procedure  can  be  summarized  in  flow  chart  Fig.  1  and  the  Algorithm 
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Fig.  1.  Flow  chart  of  APBA 


Algorithm.  In  an  imaging  experiment,  assume  we  have  performed  two  imaging  measurements. 
One  image  (radiograph)  is  the  attenuation  image  Aq  acquired  at  SID=  R\,  another  is  a  nor¬ 
malized  phase-contrast  image  I  =  M 2  ■  Id(/d) /ha  acquired  at  SID=  R \  +  Rz-  With  Aq  and  I  as 

well  as  the  initial  81,  usually  0,  as  known  inputs,  we  first  assume  /kn  =  (f/l  +  V8l^j  ■  Then 
(1).  Compute  Akn  =  \/5)  (/kn)  and  (j)  from  the  duality  transform  Eq.  (6). 
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(2).  Compute  S A  from 


SA  =  Akn  •  (1  —  P),  P  =  Aq/Akm. 


(10) 


Equations  (9)  and  (10)  are  in  fact  the  same  equations.  The  advantage  of  Eq.  (10)  over 
(9)  is  that  we  can  set  a  threshold  for  P.  We  know  P  =  Ape,coh  in  the  ideal  case  and 
Ape.coh  is  bounded  between  1  and  Aq.  The  computation  rounding  error  or  the  presence  of 
measured  noise  in  the  acquired  data  can  make  the  value  of  P  over  pass  these  bounds  in 
the  iterative  computations.  By  setting  a  threshold  upper  bound  ubd  =  1  and  lower  bound 
lbd  =  min(Ao),  the  minimum  value  ofAo,  to  P  in  the  iterative  computations,  we  can  make 
the  algorithm  more  stable.  Moreover,  if  we  know  a  better  lbd  for  Ape,coh,  other  than  the 
minimum  ofAo,  the  converging  speed  of  the  algorithm  can  be  greatly  improved. 

(3).  Compute  81  by  Fresnel  propagate  SAe11^  from  position  R\  to  R2:  Si  =  \fjx(8Ae'^\f  with 
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(4) .  Compute  /kn 


Go  to  (1 )  for  next  iteration. 


The  number  of  iterations  or  an  accuracy  measure  can  be  used  to  determine  when  to  exit  the 
program:  assuming  | -|  is  some  kind  of  norm,  that  can  effectively  reflect  the  accuracy  of  the 
retrieved  data  as  an  image,  if  ||54+i  —  5/*||  is  less  than  a  predefined  threshold  value  /j,  or  the 
iteration  step  exceeds  a  predefined  maximum  number  of  iteration  steps,  then  0  is  the  retrieved 
phase  and  the  iteration  stops;  otherwise  further  iteration  is  needed. 

An  appropriate  image  accuracy  measure  should  be  a  measure  that  can  effectively  reflect 
the  accuracy  of  the  data  AND  at  the  same  time  correctly  reflect  the  visual  perception  of  the 
data  as  an  image,  since  an  image’s  visual  perception  is  crucial  for  diagnostic  radiology.  In  the 
next  subsection  the  authors  suggest  a  measure  which  can  be  employed  as  an  appropriate  image 
accuracy  measure,  which  was  first  investigated  by  Rudin  in  [32], 


2.2.  An  Image  Accuracy  Measure 


(a)  ( b )  (c) 

Fig.  2.  Example  “Lena”  images  used  in  measuring  the  closeness  between  images,  (a)  <j)\ , 
(b)  02,  (c)  03 


A  continuous  signal  is  generally  represented  as  a  function  of  vector  variables:  /( r).  A  sam¬ 
pled  signal  will  be  represented  as  a  one-  (or  higher)  dimensional  sequence  of  real  numbers.  In 
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this  paper,  we  will  denote  the  continuous  two  dimensional  image  as  an  intensity  function  of 
two  dimensional  variables,  such  as  or  /( ?),  r  =  (x,y).  The  sampled  2-D  image  will  be 

represented  by  i  =  1 , 2,  •  •  -  ,m,  j  =  1 , 2,  -  •  •  ,n.  In  practice,  to  estimate  a  true  signal  in 

noise,  the  most  frequently  used  methods  are  based  on  the  least  squares  criteria  and  thus  an  intu¬ 
itive  measure  for  the  closeness  of  two  image  functions  /  and  g  is,  similar  to  statistical  standard 
deviation,  based  on: 


std (g,f)  :=  1 


fn(g(f)-f(f)-n)2df 


1/2 


V(0) 


(12) 


where  £2  is  the  finite  domain  of  image  functions  /  and  g,  V  (£2)  represents  the  area  of  domain 
£2,  and  fi  =  fn(g  —  f)dr/V(£2)  is  the  statistical  mean  value  of  g  —  f.  In  [32],  L.  I.  Rudin 
investigated  the  relation  of  edge  formation  of  the  2-D  digital  image  and  the  singularities  of 
the  image  function  and  pointed  out  that  the  image  intensity  function  belongs  to  the  space  of 
functions  of  bounded  total  variation.  Rudin  et  al.[33]  pointed  out  that  the  proper  norm  for 
images  is  the  total  variation  (TV)  norm,  which  is  the  L\  norm  of  the  gradient  of  the  image 
function,  and  not  the  Li  norm.  For  two  image  functions  /  and  g,  we  define  the  TV  norm  of 
g  —  f  as  the  closeness  of  the  two  images: 


T  V(g,f) 


V(£2) 


dy 


V(£2) 


(13) 


where  V  is  the  gradient  operator.  Since  std  is  the  form  of  Li  norm,  it  is  not  suitable  as  a 
measure  to  represent  the  closeness  between  two  images.  For  example,  for  the  three  “Lena” 
image  functions  shown  in  Fig.  2,  the  std  measures  of  (f)\  —  02  and  01  —  03  are  std(0i ,  0?)  =  0.206, 
std(0i ,  03)  =  0.472  respectively.  But  the  image  0 3  is  much  more  closer  to  0 1  than  0 2  does  in 
visual  perception.  This  is  because  the  digital  representation  of  an  image  depends  not  only  on 
the  pixel  values,  it  also  depends  on  the  contrast  changes  between  neighbor  pixels.  An  image’s 
visual  perception  is  crucial  for  diagnostic  radiology.  This  contrast  changes  between  neighbor 
pixels  can  better  be  represented  by  the  gradient  changes  of  the  image  functions.  For  example, 
the  TV  measures  of  0i  —  0 2  and  0 1  —  0 3  are  TV(0i,02)  =  0.0247  and  TV(0i,03)  =  0.0138 
respectively,  more  appropriate  in  reflecting  the  visual  perception.  In  this  paper,  we  will  use  the 
TV  norm  (13)  as  the  measure  for  closeness  between  two  compared  image  functions. 

Note  that  the  TV  norm  between  two  image  functions,  say  g  and  /,  equals  0  if  and  only 
if  g  differs  from  /  by  a  constant.  This  feature  does  not  affect  its  appropriateness  for  phase 
retrieval  since  it  is  well  known  that  the  recovered  phase  0  is  unique  up  to  a  constant  with  given 
information  about  the  attenuation-map  ,4q  and  the  phase-contrast  intensity-map  I. 


3.  Simulation  Tests 

In  order  to  investigate  the  performance  of  the  algorithm  constructed  above,  we  perform  com¬ 
puter  simulations  in  this  section.  In  Section  3.1,  we  first  construct  a  breast  tissue  model  that 
represents  the  phase-shifts  and  attenuation  of  breast  tissues  and  embedded  microcalcifications 
for  different  x-ray  energies.  In  our  simulation  tests,  the  distances  of  source  point  to  object,  Ri, 
and  object  to  detector,  R2,  are  set  to  26  inches  (0.66  m)  respectively.  In  this  way  the  magni¬ 
fication  factor  M  =  (R\  +  Ry) / R\  is  equal  to  2.  For  convenience,  the  incident  x-ray  intensity 
/in  at  R\  is  set  to  M 2  (4).  We  compare  the  performance  of  our  algorithm,  APBA,  with  two 
other  widely  used  algorithms:  the  Gerchberg-Saxton  (GS)  algorithm  (Section  3.2)  and  the  TIE 
algorithm  (Section  3.3). 
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3.1.  A  breast  specimen  model 

The  tissue’s  phase-shift  and  attenuation  are  determined  not  only  by  the  tissue’s  physical  com¬ 
position,  but  also  by  the  x-ray  energy.  With  different  x-ray  energy,  the  same  tissue  has  different 
phase-shift  and  attenuation  change.  Simulation  models  used  in  literature  often  do  not  incorpo¬ 
rate  these  changes.  In  this  subsection,  we  construct  a  breast  specimen  model  that  can  represent 
tissue  attenuation  and  phase  shifts  according  to  employed  x-ray  energies  as  well  as  tissue’s 
compositions. 


(a) 


(b) 


(c) 


Fig.  3.  Image  manifest  of  (a)  A^e  coh,  (b)  Aj^  and  (c)  Ag  when  x-ray  energy  equals  35.5  keV. 


(°)  (b) 


(c) 


Fig.  4.  Profiles  of  A^e  coh,  the  solid  lines,  and  Aj^,  the  dotted  lines,  when  x-ray  energy 
equals  (a)  18.5  keV,  (b)  35.5  keV  and  (c)  59.5  keV. 


In  our  model  the  tissue  has  two  physical  compositions:  the  50%  glandular-50%  adipose 
breast  tissue  and  the  microcalcifications.  In  order  to  simulate  the  morphological  aspects  of 
breast  tissues,  we  adopted  a  radiograph  of  a  breast  lumpectomy  specimen  with  pixel  values 
rescaled  and  the  metal  localization  wire  removed  by  replacing  the  pixel  value  at  wire  position 
with  a  mean  pixel  value  at  near  by  positions.  It  is  difficult  to  remove  all  the  residual  trace  of 
the  wire  this  way  as  can  be  seen  from  the  following  image  display  especially  for  the  phase 
contrast  image  (Fig.  5(c)).  In  the  phase  contrast  image  (Fig.  5(c)),  the  small  residual  variation 
from  the  original  wire-track  really  got  enhanced.  The  linear  attenuation  coefficients  and  elec¬ 
tron  densities  for  the  50%  glandular-50%  adipose  breast  tissues  are  computed  from  the  tissue’s 
elemental  composition  and  the  interpolated  elemental  mass  attenuation  coefficients  from  the 
standard  reference  in  the  mammographic  radiation  dosimetry  [34,  35],  Moreover,  each  mass 
attenuation  coefficient  is  further  broken  down  to  two  components:  one  from  x-ray  photoelec¬ 
tric  absorption  and  coherent  scattering,  and  another  from  incoherent  scattering.  In  this  way, 
the  total  attenuation  is  partitioned  as  a  product  of  A“c  coh  and  as  defined  in  Eq.  (8)  above. 


#126430  -  $15.00  USD  Received  2  Apr  2010;  revised  1 1  Jun  2010;  accepted  9  Jul  2010;  published  15  Jul  2010 
(C)  2010  OS  A  19  July  2010/ Vol.  18,  No.  15  /  OPTICS  EXPRESS  16082 


In  addition,  to  simulate  the  microcalcifications  in  breast,  four  small  ellipsoids  of  calcium  car¬ 
bonate  (CaC03)  are  embedded  in  the  specimen  model.  The  diameter  of  the  ellipsoids  can  be 
adjusted  in  simulating  different  size  of  the  microcalcifications.  In  the  following  simulations,  to 
test  phase-contrast  sensitivity,  we  set  the  diameters  of  the  four  ellipsoids  to  10,  5,  10,  and  5 
microns  in  x-ray  direction  and  300,  200,  300,  200  microns  in  detector  plane  respectively. 

The  attenuation  image  A 5  of  the  specimen  model  simulated  with  35.5  keV  x-ray,  and  its  two 
corresponding  partition  images  Aj^  and  Apecoh,  are  shown  in  Fig.  3.  For  a  comparison,  the 
profiles  of  A^  and  A^e  coh  along  the  line  passing  through  the  microcalcifications  are  shown  in 
Fig.  4  simulated  with  x-ray  energy  equals  18.5,  35.5  and  59.5  keV,  respectively.  We  can  see  that 
the  contribution  of  ApC  coh  to  the  total  attenuation  A5  gets  smaller  when  x-ray  energy  is  getting 
higher.  Especially,  when  x-ray  energy  equals  59.5  keV,  the  contribution  of  A^e  coh  for  the  soft 
tissue  can  almost  be  neglect,  as  is  expected. 

3.2.  Comparison  with  the  GS  Algorithm 


(a)  (6)  (c) 


Fig.  5.  Image  representation  of  the  inputs  generated  from  the  simulation  model  and  Fresnel 
propagation,  (a)  the  phase  map  <j>;  (b)  the  attenuation  map  Aq;  and  (c)  the  normalized  Fres¬ 
nel  propagated  phase  contrast  image  I  with  object  to  detector  distance  JO  =  26  in  (0.66  m). 

The  GS  algorithm  is  an  iterative  algorithm  for  phase  retrievals  from  a  pair  of  images  at 
two  planes  related  by  the  Fourier  transform.  For  details  readers  are  referred  to  [26],  The  GS 
algorithm  is  a  classical  algorithm  which  is  widely  used  in  electron  microscopy,  wave  front 
sensing,  astronomy,  crystallography,  and  many  other  fields  involving  phase  recovery  [27,  28, 
36], 

By  replacing  the  Fourier  Transform  in  the  GS  algorithm  with  the  Fresnel  transform  of 
Eq.  (11),  we  get  a  modified  version  of  the  GS  algorithm  extended  to  the  Fresnel  diffraction 
regime.  Our  previous  simulation  tests  showed  that  this  modified  GS  algorithm  converges  very 
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slow  for  object-detector  distance  Rs  ~  1  m,  and  converges  faster  for  larger  /G,  such  as  images 
acquired  at  synchrotron  beam  lines.  But  this  is  generally  not  applicable  to  the  held  of  clinical 
imaging,  due  to  the  physical  constraints  such  as  compact  sizes  of  hospital  rooms. 

In  our  simulation  tests,  we  compare  the  performance  of  APBA  and  that  of  the  GS  algorithm. 
The  photon  energy  of  the  point  x-ray  source  is  set  to  35.5  keV,  and  the  distances  from  the 
source  to  object  and  object  to  detector  are  set  to  R\  =  R2  =  26  inches  (0.66  m)  respectively. 
For  convenience,  the  incident  x-ray  intensity  Tm  at  R\  is  set  to  M2,  where  M  =  (R\  +  Rf)/R\ 
is  the  magnification  factor.  The  phase  map  (j)  and  the  attenuation  A5  are  generated  from  our 
phantom  model  for  35.5  keV  x-ray.  Fig.  5  shows  the  simulated  images  of  the  phase  map  (j).  the 
attenuation  image  A5  and  the  phase-contrast  image  I. 


(a)  (b)  (c) 

Fig.  6.  Comparison  of  the  performance  of  the  GS  algorithm  and  APBA.  (a)  plot  of  the 
accuracy  measures  with  respect  to  iteration  steps.  The  plot  with  solid  line  represents  the 
APBA.  The  one  with  dashed  line  represents  the  GS  algorithm;  (b)  recovered  phase  map 
with  the  GS  algorithm  after  100  iterations;  (c)  recovered  phase  map  with  APBA  after  100 
iterations. 

In  the  simulation  test,  the  iteration  for  our  attenuation-partition  based  algorithm  (APBA) 
and  the  GS  algorithm  are  performed  100  steps.  The  corresponding  recovered  phase  images  are 
shown  in  Fig.  6(b)  and  (c).  In  Fig.  6(a),  the  solid  line  represents  the  change  of  total  variation 
(TV)  of  the  retrieved  phase  using  attenuation-partition  based  algorithm  with  respect  to  the  it¬ 
eration  steps.  The  dashed  line  represents  the  change  of  TV  of  the  retrieved  phase  using  the 
(modified)  GS  algorithm  with  respect  to  iteration  steps.  We  can  see  that  the  converging  speed 
of  the  (modified)  GS  algorithm  is  much  slower  than  that  of  attenuation-partition  based  algo¬ 
rithm  (the  change  of  TV  of  APBA  from  step  1  to  step  100  is  1 .04 E  —  3,  almost  33  times  greater 
than  that,  3 . 172?  —  5,  of  the  GS  algorithm.)  In  addition,  from  the  visual  perception  point  of  view, 
we  see  that  the  phase  map  retrieved  with  the  attenuation-partition  based  algorithm  (APBA)  is 
much  better  than  that  retrieved  with  the  (modified)  GS  algorithm. 

The  main  difference  between  APBA  and  the  GS  algorithm  is  that  APBA  takes  the  advan¬ 
tage  of  the  phase-attenuation  correlation  property,  although  the  extent  of  this  correlation  is  not 
known  in  priori ,  but  the  GS  algorithm  does  not.  From  this  example  we  see  that  the  phase- 
attenuation  correlation  is  a  very  important  information  that  shouldn’t  be  neglected  in  the  algo¬ 
rithm  development  for  phase  retrieval. 

3.3.  Comparison  with  the  Transport  of  Intensity  (TIE)  algorithm 

We  have  mentioned  the  transport  of  intensity  equation  in  Eq.  (3)  in  Section  1.  Since  Teague 
proposed  the  TIE  algorithm  for  phase  retrieval  in  1983  [21],  numerous  phase  retrieval  algo¬ 
rithms  have  been  suggested  on  how  to  effectively  search  the  numerical  solution  of  the  TIE 
[4,  20,  25,  37,  38,  39,  40].  Among  the  methods  of  solving  the  above  TIE  for  the  phase  map, 
the  one  based  on  Fast  Fourier  Transform,  proposed  by  Nugent  et  al.  [4],  and  Allen  and  Oxley 
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in  [25],  is  the  most  widely  used.  In  the  form  of  pseudo-differential  operator,  the  solution  phase 
map  0  is  given  by 


<Hr) 


v[v-2(/-a2)] 


^5 


(14) 


for  the  given  normalized  phase-contrast  image  I  and  attenuation  image  A5.  Here  V  is  the  gradi¬ 
ent  operator,  and  V  2  is  the  inverse  Laplacian  operator. 

The  advantage  of  this  algorithm  is  that  it  does  not  require  the  boundary  information  in  solv¬ 
ing  the  TIE  (assuming  the  image  data  is  periodic);  it  is  a  deterministic  method  and  thus  the 
algorithm  is  fast  and  accurate  comparing  to  most  iterative  algorithms.  In  this  section  we  com¬ 
pare  the  performance  of  the  TIE  algorithm  and  that  of  APBA  for  two  kind  of  cases:  first  for  the 
ideal  case  without  any  noise  and  any  image  misalignment,  then  for  cases  simulating  practical 
situation  with  x-ray  imaging  noise  and  possible  image  misalignment.  In  these  simulation  tests, 
the  imaging  geometries  are  the  same  as  in  the  previous  subsection,  and  x-ray  energy  is  again 
35.5  keV  . 

For  the  ideal  case  without  any  noise  and  any  image  misalignment,  the  performance  compar¬ 
ison  results  are  shown  in  Fig.  7.  For  the  ideal  case  the  TIE  algorithm  is  accurate  both  in  TV 
measure  and  visual  perception.  APBA  can  also  achieves  this  accuracy  but  needs  1110  iteration 
steps  to  have  its  TV  measure  of  0.00215226,  an  error  smaller  than  that  of  0.00215269  with  the 
TIE  algorithm. 

However,  the  real  test  lies  in  the  performance  for  the  cases  simulating  practical  situation 
with  x-ray  imaging  noise  and  possible  image  misalignment.  Obviously  only  the  performance  in 
these  cases  really  matters  in  phase  contrast  imaging  applications,  especially  for  clinical  imaging 
applications  where  the  imposed  radiation  limits  dictates  existence  of  substantial  x-ray  quan¬ 
tum  noise  in  acquired  images.  Implemented  in  the  Fourier  space,  the  inverse  Laplacian  V  2 
in  Eq.  (14)  is  singular  at  zero  spatial  frequency.  This  singularity  will  amplify  the  noise  in  the 
images  and  result  in  failure  of  accurate  phase  retrieval  for  the  TIE  algorithm.  To  overcome  this 
difficulty,  some  kind  of  regularization  must  be  used.  The  most  widely  used  regularization  is 
Tikhonov  regularization,  which  replaces  V~2  by  |u|2/(|u|4  +  ff2),  for  some  “favorable”  con¬ 
stant  parameter  K,  called  the  Tikhonov  regularization  parameter.  In  this  regularization  scheme, 
the  singularity  is  regularized,  and  the  favorable  parameter  K  means  the  retrieved  phase  (j)K  corre¬ 
sponding  this  K  is  as  close  to  the  true  phase  0  as  possible.  Roughly  speaking,  the  regularization 


Fig.  7.  Comparison  of  APBA  and  the  TIE  algorithm  with  pure  data,  (a)  plot  of  the  accuracy 
measure  with  respect  to  iteration  steps.  The  plot  with  solid  line  represents  the  APBA.  The 
one  with  dashed  line  represents  the  TIE  algorithm.  It  needs  1110  steps  for  the  TV  measure, 
0.00215226,  of  APBA  to  achieve  to  the  TV  measure,  0.00215269,  of  TIE;  (b)  recovered 
phase  using  the  TIE  algorithm;  (c)  recovered  phase  using  APBA  after  1500  iteration  steps. 
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parameter  K  is  inversely  proportional  to  the  images  signal-noise  ratio.  Two  problems,  however, 
arise  to  this  regularization.  First,  the  true  phase  <j>  is  not  known  in  real  situations.  So  the  regu¬ 
larization  parameter  K  is  not  a  priori,  which  makes  it  difficult  in  practical  applications.  Second, 
this  Tikhonov  regularization  is  based  on  finding  a  stable  solution  to  Ax  =  y,  in  Hilbert  spaces 
X,  Y,  by  solving  the  minimization  problem 

xK  =  argmin  ||Ax  — y||y  +  k:||x|||  .  (15) 

xex 

It  is  Li  norm  dependent.  Since  the  proper  norm  for  image  data  is  the  total  variation  (TV)  norm 
[33],  a  favorable  solution  under  the  Tikhonov  regularization  principle  can  not  be  guaranteed 
a  best  solution  in  visual  perception.  Moreover,  for  relatively  noisy  acquired  images,  the  TIE- 
based  algorithm,  even  with  Tikhonov  regularization,  often  failed  to  retrieve  the  phase  maps[l]. 


(<*)  (e)  (/) 


Fig.  8.  Comparison  of  the  TIE  algorithm  and  APBA  with  noise  added,  (a)  True  phase  map 
</);  (b)  attenuation  map  A?;  (c)  the  normalized  phase  contrast  image  /;  (d)  recovered  phase 
map  with  the  TIE  algorithm,  no  Tikhonov  regularization  is  used;  (e)  recovered  phase  map 
with  the  TIE  algorithm  with  Tikhonov  regularization;  (f)  recovered  phase  map  with  APBA 
after  10  iteration  steps.  In  the  simulation,  the  acquired  data  is  assumed  to  have  a  level  of 
8i,  =  0.03%  detector  noise  and  one  pixel  misalignment  between  A^  and  /  horizontally. 

In  the  following,  we  will  compare  the  performance  of  APBA  and  the  TIE-algorithm  when  noise 
is  present.  In  the  practice  of  phase  retrieval,  there  are  generally  two  kinds  of  image  data  errors. 
One  is  the  noise  associated  with  image  acquisitions,  including  the  quantum  noise  of  x-ray  pho¬ 
ton  detections  and  detector  electron  noise.  We  assume  the  quantum  noise  dominates  as  is  the 
case  in  most  imaging  applications.  The  other  is  the  error  caused  by  the  misalignment  between 
the  attenuation  map  Aq  and  phase-contrast  image  I.  This  is  because  usually  the  attenuation 
image  and  corresponding  phase  contrast  image  are  generally  acquired  in  two  separate  x-ray 
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exposures.  The  misalignment  could  be  resulted  from  the  shift  or  tilting  of  the  object  or  detector 
between  the  exposures.  In  the  simulation,  we  associate  each  pixel  value  of  an  image  a  photon 
count  V,  so  that  P(i.j)  =c-N(i,j),  where  c  is  a  constant.  Assuming  the  noise  has  a  Poisson  dis¬ 
tribution,  with  variance  a2  =  N  at  each  pixel.  In  the  simulation  we  assign  a  background  noise 
level  for  each  simulated  image.  This  background  noise  level  is  defined  as  the  ratio  8b  =  (?b/Nb 
corresponding  to  the  direct  exposure  area  (where  A2  =  1)  outside  the  object  in  background.  The 
Poisson  statistics  dictates  that  Nb=  1  /§£  and  in  this  way  the  photon  count  N (i.  j)  can  be  deter¬ 
mined  accordingly  at  each  pixel.  Once  N(i,j)  is  determined,  the  statistical  errors  at  each  pixel 
can  be  assigned  using  a  computer  simulated  random  Poisson  distribution  generator  with  mean 
corresponding  to  the  photon  counts  N(i,j).  In  the  simulations  below,  the  background  noise 
level  8b  is  set  to  0.0003.  The  images  with  the  noise  added  are  shown  in  Fig.  8(b)  and  (c).  The 
quality  of  an  image  depends  not  only  on  the  noise  level  but  also  on  the  extent  of  image  contrast 
change.  With  the  assumption  P(i,j)  =  c -N(i,j)  above,  one  can  easily  see  that  8/,  =  (Jpb/Pb  is 
the  statistical  coefficient  of  variation  in  absorbed  photon  numbers  in  background.  8()i  =  OpJPoi, 
the  coefficient  of  variation  of  the  object  image,  on  the  other  hand,  is  the  structural  coefficient  of 
variation  of  the  sampled  image,  which  reflects  the  image’s  normalized  extent  of  image  contrast 
change.  We  will  use  the  ratio  8ol  /  8b  to  reflect  the  image  quality.  The  larger  the  ratio,  the  higher 
the  image  quality.  In  our  simulation  models,  when  8b  =  0.0003,  the  corresponding  ratio  8t)i/8b 
for  attenuation  A5  and  phase-contrast  image  /,  Fig.  8(b)  and  (c),  are  1.67  and  6.35  respectively. 
Because  of  the  phase-contrast  effect,  the  image  quality  of  the  phase-contrast  image  I  is  higher 
than  the  attenuation  image  A5  although  they  have  the  same  background  noise  level. 

Three  simulation  tests  are  performed  in  the  comparison:  easel:  assume  the  acquired  data  Aq 
and  I  has  noise  present  but  perfectly  aligned;  case  2:  assume  acquired  data  has  no  noise  but 
has  one  pixel  misalignment  horizontally;  case  3:  assume  acquired  data  has  combined  detector 
noise  as  well  as  one  pixel  misalignment.  One  bias  in  simulation  for  the  TIE-algorithm  should 
be  mentioned:  with  the  known  true  phase  value,  a  favorable  Tikhonov  regularization  parameter 
K  can  be  searched,  but  in  practice  this  search  is  hardly  feasible.  In  each  case  mentioned  above, 
three  phase  retrievals  are  performed:  1 .  using  the  TIE-algorithm  without  Tikhonov  regulariza¬ 
tion;  2.  using  the  TIE-algorithm  and  the  favorable  Tikhonov  regularization  parameter;  3.  using 
APBA  with  10  step  iteration.  The  Total  variation  (TV)  of  the  results  are  listed  in  Tab.  1  and 
the  recovered  phase  images  for  case  3  are  shown  in  Fig.  8(d)  -  (f).  The  results  using  Tikhonov 
regularization  are  better  than  those  not  using  Tikhonov  regularization  but  worse  than  those  us¬ 
ing  APBA.  The  influence  of  misalignment  to  APBA  is  little  but  disaster  to  the  TIE  algorithm. 
From  the  profile,  shown  in  Fig.  9,  along  a  line  passing  through  the  microcalcifications  we  can 
see  that  the  values  of  the  phase  recovered  from  APBA  is  very  close  to  the  true  phase  value,  but 
the  values  of  the  phase  recovered  from  the  TIE  algorithm  are  distorted. 

4.  Discussion  and  Conclusion 


Table  1 .  TV  comparison  of  the  TIE  algorithm  and  APBA.  In  the  table,  K  is  the  Tikhonov 
regularization  parameter,  A  represents  the  sampling  step-size  in  FT-space. 


TV  = 

Case  1 

Case  2 

Case  3 

TIE:  no  Tikhonov 

0.0406 

0.0808 

0.0886 

TIE:  with  Tikhonov 

K  =  0.7828A2 

0.0399 

K  =  4.4421 A2 

0.153 

K  =  4.1076A2 

0.0379 

APBA 

0.0283 

0.0086 

0.0290 
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Fig.  9.  Profiles,  along  a  line  passing  through  the  microcalcifications,  of  the  recovered  phase 
using  APBA,  the  solid  line,  and  using  the  TIE  algorithm  with  Tikhonov  regularization,  the 
dashed  line,  in  Case  3.  The  dash-dotted  line  is  the  true  phase. 


In  above  analysis,  we  assumed  that  the  x-ray  source  is  a  quasi-monochromatic  point  source.  In 
practice,  one  often  employs  conventional  incoherent  and  polychromatic  sources  such  as  x-ray 
tubes  for  imaging.  In  the  experiments  with  an  x-ray  tube  source,  the  previous  formula  Eqs.  (2- 
3)  of  the  phase  contrast  image  intensities  should  be  modified.  Since  in  the  APBA  method  the 
duality  transform  is  derived  based  on  Eq.  (3),  hence  the  Eq.  (6)  should  be  modified  accordingly. 
In  our  previous  works  we  have  studied  this  problem  [20].  With  the  Wigner  function  based  phase 
space  formalism,  we  have  proved  that  the  coherence  degree  of  a  finite-size  focal  spot  can  be 
accounted  for  by  the  optical  transfer  function  OTFq  jj  (h/M)  for  the  geometrical  unsharpness 
associated  with  the  finite-size  source  [20]: 


otfg. 
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where  /spot(4)  is  the  intensity  distribution  of  the  focal  spot.  We  found  the  generalized  TIE 
equation  with  an  x-ray  tube  source,  that  is,  the  x-ray  intensity  at  the  detector  is  given  by  [20]: 
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where  operator  (g>  denotes  the  convolution,  A 2  is  the  total  attenuation  of  the  imaged  object,  0  is 
the  spectrally  averaged  phase-shift  caused  by  the  object,  and  ( •  )  means  the  spectral  average. 
Compare  above  equation  with  Eq.  (3)  and  it  is  clear  that  the  TIE-based  phase  retrieval  method 
needs  only  two  modifications:  (i).  Fourier  space  de-convolution  of  the  measured  intensity  from 
OTFq  jj  (u/M),  (ii).  Replacing  wavelength  with  the  spectral-averaged  (A2)/(A).  In  the  same 
fashion,  the  duality  transform  D  defined  in  Eq.  (6)  should  be  modified  with  (i).  Fourier  space  de- 
convolution  of  the  measured  intensity  from  OTFqjj.{u/M)\  (ii).  A  replacement  of  the  parameter 
k  defined  in  Eq.(7)  with  the  spectral-average 
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Otherwise,  the  APBA  flow  chart  is  the  same  as  that  for  the  case  with  a  quasi-monochromatic 
point  source. 

Phase  retrieval  is  a  crucial  step  for  quantitative  imaging  such  as  reconstructing  the  3-D  dis¬ 
tribution  of  tissue  linear  attenuation  coefficients  and  refraction  indices.  However,  images  ac¬ 
quired  in  medical  applications  are  relatively  noisy,  due  to  the  radiation  dose  constraints,  with 
low  phase  contrast  effect,  due  to  physical  constraints  such  as  compact  sizes  of  hospital  rooms. 
An  phase  retrieval  algorithm  which  is  robust  to  noise  is  necessary  for  potential  medical  phase 
contrast  imaging.  In  [1],  the  authors  developed  an  algorithm,  called  attenuation-partition  based 
phase  retrieval  algorithm.  It  is  an  iterative  algorithm  which  takes  advantage  of  the  correlation 
between  the  attenuation  and  phase-shift.  The  phase  retrieval  results  from  experimental  images 
show  that  this  algorithm  is  fast  and  robust  [1],  In  this  work,  we  systematically  compared  the 
performance  of  this  algorithm  with  other  two  widely  used  phase  retrieval  algorithms,  namely 
the  Gerchberg-Saxton  (GS)  algorithm  and  the  Transport  of  Intensity  Equation  (TIE)  algorithm. 
The  systematical  comparison  is  conducted  by  analyzing  phase  retrieval  performances  with  a 
digital  breast  specimen  model.  We  show  that  the  proposed  algorithm  converges  faster  than  the 
GS  algorithm  in  the  Fresnel  diffraction  regime,  and  is  more  robust  against  image  noise  than 
the  TIE  algorithm.  These  results  suggest  the  significance  of  the  proposed  algorithm  for  future 
medical  applications  with  the  x-ray  phase  contrast  imaging  technique. 
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In  the  inline  phase  contrast  X-ray  tomography  the  reconstructed  apparent  linear  attenuation  coefficient 
values  may  be  greatly  larger  than  sample’s  linear  attenuation  coefficients  or  even  be  negative.  In  this 
work  we  present  a  general  formula  to  quantitatively  relate  the  apparent  linear  attenuation  coefficient 
values  in  cone-beam  phase  contrast  tomography  to  sample’s  linear  attenuation  coefficients  and  refractive 
indices.  This  formula  overcomes  the  gross  inaccuracy  of  the  existing  formula  in  the  literature  in  analyzing 
high-resolution  phase  contrast  tomography,  and  it  will  be  useful  for  correctly  interpreting  and  quantify¬ 
ing  the  apparent  linear  attenuation  coefficients  in  cone-beam  X-ray  phase  contrast  tomography. 
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For  many  applications  the  X-ray  attenuation  contrast,  such  as 
mass-density  differences  between  different  material  components, 
is  too  small  to  generate  enough  image  contrast  with  conventional 
X-ray  tomography.  Taking  advantage  of  highly  sensitive  X-ray 
phase  shifts  generated  from  X-ray  coherent  scattering,  phase  con¬ 
trast  X-ray  tomography  finds  broad  applications  in  material  sci¬ 
ence  and  engineering,  biomedical  imaging  and  many  other 
scientific  fields  [1-3].  In  phase  contrast  X-ray  tomography  one  first 
makes  phase  contrast  angular  projections,  in  which  the  refraction 
and  diffraction  of  phase  shifted  X-rays  form  the  edge-enhancing 
interference  fringes  at  interfaces  of  different  material  components. 
From  these  angular  projections  3-D  tomograms  are  then  recon¬ 
structed  using  the  standard  tomography  reconstruction  methods 
such  as  the  filtered  backprojection  methods.  The  reconstructed 
phase-sensitive  tomograms  exhibit  enhanced  interfaces  and  layers 
of  different  materials  and  tissues  [1-3].  However,  these  tomograms 
are  not  the  maps  of  the  linear  attenuation  coefficients  (LACs)  of  the 
sample,  rather  they  are  the  maps  of  the  reconstructed  apparent 
LACs,  which  may  present  anomalously  large  or  even  negative 
apparent  LAC  values  at  interfaces  between  different  materials 
and  tissues.  These  “artifacts”  may  cause  faulty  interpretation  of 
sample  structures,  and  impede  even  qualitative  characterizations 
for  tissues  and  materials  [1,2, 4-6].  As  a  result,  for  example,  the 
apparent  LAC  values  of  small  bronchi  were  found  lower  than  that 
of  air  in  phase  contrast  tomography  experiments  [2],  In  order  to 
correctly  interpret  and  quantify  the  mixed  contrast  exhibited  in 
phase  contrast  tomography,  it  is  necessary  to  find  out  quantitative 
relations  between  the  apparent  LAC,  on  the  one  hand,  and  the  ac¬ 
tual  LAC  and  refractive  index,  on  the  other.  The  previous  studies 
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found  that  the  apparent  LAC  /iapparent(ro)  in  X-ray  phase  contrast 
tomography  is  given  by  [7-9]: 

/W?o)  =  H(fo)-g-  (!  +  £  +  !)Wo).  (1) 

In  this  equation  the  first  additive  term  /t(r0)denotes  the  sample’s 
actual  LAC  value  at  position  r0,  and  the  second  additive  term  is  pro¬ 
portional  to  the  three-dimensional  Laplacian  of  the  sample’s  refrac¬ 
tive  index  decrement  Arefr(r0)  [7-9].  In  Eq.  (1)  R2  denotes  the 
sample-detector  distance  and  M  the  magnification  factor  employed 
in  angular  projections.  In  the  literature  this  second  term  is  used  to 
explain  the  edge  enhancement  and  the  anomalous  values  of  large 
apparent  LAC  [7-9]. 

However,  in  high-resolution  phase  contrast  tomography  with 
micrometer-scale  resolutions  the  applicability  of  Eq.  (1)  is  called 
into  question.  Researchers  recently  found  that  Eq.  (1)  is  invalid  in 
some  synchrotron  radiation  microtomography  experiments  [10]. 
In  fact,  in  the  derivation  of  Eq.  (1)  the  phase-sensitive  projections 
were  modeled  based  on  the  so-called  Transport  of  Intensity  Equa¬ 
tion  (TIE)  [7-9].  The  TIE  describes  how  the  phase-shifted  X-rays  re¬ 
fract  and  diffract  to  form  interference  fringes  at  boundaries  of 
material  components  [11,12].  However,  we  noted  that  TIE  can  only 
describe  the  moderate-resolution  phase-sensitive  projections  with 
nAR2/(4Ma2)  «  1,  where  a  denotes  the  finest  size  resolved.  In 
cases  of  high-resolution  projections  with  micrometer-scale  resolu¬ 
tions,  the  fine  size  resolved  could  be  so  small  such  that 
nAR2/4Mo2  ~  1  or  larger,  then  TIE  becomes  inapplicable  and  Eq. 
(1)  cannot  quantify  the  apparent  LAC  values  observed  in  high- 
resolution  phase  contrast  tomography. 

To  overcome  the  limitations  of  Eq.  (1),  we  derived  a  novel  and 
general  formula  of  apparent  LAC  values  in  phase  contrast  X-ray 
tomography.  Instead  of  using  the  TIE,  in  the  derivation  we 
employed  the  general  phase-sensitive  projection  formula  that  is 
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applicable  for  high-resolution  projections  as  well  [13-14].  We  then 
applied  the  standard  Feldkamp-Davis-Kresss  (FDK)  cone-beam 
reconstruction  algorithm  to  the  ray-sums  of  phase-sensitive  angu¬ 
lar  projections  [15].  We  found  that  the  apparent  LAC  value 
H, pparent(ro)  is  given  by 

.  (  aR2  _2  \  ,  47t 

ftpparent(ro)  =  COS  3D  j  fi(r0)  -  — 

xsin(iSiv'D)Arcfr(?o)’  (2) 

where  the  three-dimensional  Laplacian  operator  \l\D  =  {d2 1  dx\+ 
d2 1  dx\  +  d2  /  dx\) .  The  operators  sin  (^(/iR2/4nM)'VlDsJ  and  cos((xR2 
/4nM)\7lD)  in  above  equation  are  the  pseudo-differential  operators. 
The  action  of  a  pseudo-differential  operator  such  as  sin(cV3D)  on  a 
function  g(r0)  is  defined  as  sin(cV3D)g(r0)  =  f  f  exp  (i2n 
(?  —  r0)  ■  /)  ■  sin  f-4n2cf2'jg(f)d3rd3f.  Note  that  Eq.  (2)  is  indeed  a 

generalization  of  Eq.  (1).  For  experiments  with  moderate-resolu¬ 
tions  or  short  sample-detector  distances  such  that 
nlR2/4Ma2  «  1,  we  then  have  cos((XR2/4nM)VlD)  rs  1  and 
sin((AR2/47tM)V3D)  ~  (/,R2/4nM)V2n.  Thereby  the  first  term  on 
the  RHS  of  Eq.  (2)  reduces  to  the  actual  LAC,  the  second  term  to  a 
scaled  Laplacian  of  Arefr(r0)  and  Eq.  (2)  recovers  Eq.  (1),  as  is 
expected. 

In  order  to  validate  Eq.  (2),  we  performed  numerical  simula¬ 
tions.  In  the  simulations  the  X-ray  source  was  a  point  source  of 
60  keV  X-rays.  The  LACs  and  refractive  indices  of  the  numerical 
phantom  were  constructed  as  the  sums  of  the  LACs  and  refractive 
indices  of  six  spheres  of  light  elements  with  different  radii  and  cen¬ 
ter  positions.  For  each  sphere  the  mass  densities  were  designed  to 
decrease  radially  and  diminish  rapidly  at  the  sphere’s  boundary  for 
generating  phase  contrast.  Specifically  we  set  /<(r0)  =  0.215  x 
El,/n(?o)(l/cm)  and  Arefr(r0)  =  (7.548  x  10  s)  x  ELM^o), 
where  /„(r0)  =  ^R2n-(r0-  fnf/4Rn  if  f(r0-r„)|^R„  and 
fn(/o)  =  0  otherwise.  Here  r„  and  R„  (n  =  1,2,.. .,  6)  denote  the  cen¬ 
ter  position  and  radius  of  the  respective  sphere.  In  this  way  the 
LACs  and  refractive  indices  change  continuously  without  any  steps 


Fig.  1.  Reconstructed  apparent  LAC  map  of  a  phantom  slice  from  the  phase- 
sensitive  cone-beam  tomography  with  mixed  contrast. 


but  vary  rapidly  at  the  sphere  boundaries.  Hence  strong  phase  con¬ 
trast  would  manifest  around  the  boundaries  of  the  spheres.  The 
detector  has  512  x  512  pixels  with  a  sampling  pitch  of  2  pm.  In 
the  simulated  scan  the  source  traced  a  circle  orbit,  and  the  source 
and  detector  rotated  together  over  an  angular  range  of  360°  in  0.5° 
steps.  In  order  to  simulate  X-ray  Fresnel  diffractions  involved  in 
the  angular  projections,  one  needs  to  compute  X-ray  attenuation 
values  and  phase  shifts  of  the  exit  rays  in  each  of  the  projections. 
In  general  one  needs  to  use  the  multislice  propagation  method  to 
account  for  X-ray  Fresnel  diffractions  inside  the  sample  [16].  How¬ 
ever,  the  settings  described  above  satisfy  the  so-called  thin  object 
condition  T  <  p2 / 2a0,  where  T  denotes  the  sample’s  thickness,  p  is 
the  resolution  and  X0  the  wavelength  of  60  keV  X-ray  [16].  This 
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Fig.  2.  Profiles  of  the  apparent  LACs  along  the  vertical  line  passing  the  center  of  the  central  sphere,  (a)  Profile  measured  from  Fig.  1.  (b)  Profile  calculated  based  on  Eq.  (1). 
(c)  Profile  calculated  based  on  the  general  formula  Eq.  (2). 
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Table  1 

A  comparison  of  the  fringe  peak-peak  differences  in  Fig.  2(a)-(c). 


Fringe  set 

Fringe  set 

Fringe  set 

Fringe  set 

Fringe  set 

Fringe  set 

Fringe  set 

Fringe  set 

1 

2 

3 

4 

5 

6 

7 

8 

Peak-peak  difference  of  apparent  LACs  in  Fig.  2(a)  (1/cm) 

4.648 

4.979 

6.098 

8.860 

8.860 

6.098 

4.979 

4.648 

Error  of  calculated  peak-peak  difference  based  on  Eq.  (1) 

429 

346 

199 

136 

136 

199 

346 

429 

(Fig.  2(b))  (%) 

Error  of  calculated  peak-peak  difference  based  on  Eq.  (2) 

15.0 

17.5 

13.0 

6.9 

6.9 

13.0 

17.5 

15.0 

(Fig.  2(c))  (%) 


condition  allows  us  to  compute  the  X-ray  attenuation  values  and 
phase  shifts  of  the  exit  X-rays  as  the  ray  integrals  of  the  sample’s 
LACs  and  refractive  index  decrements.  In  this  simulation  a  ray¬ 
tracing  algorithm  was  used  to  compute  the  ray  integrals  such  as 
/ /.idl  and  /  Arefrdl  in  each  of  the  angular  projections.  The  exit 
X-ray  attenuation  and  phase  shift  were  determined  as  A 2  = 
exp  (—  f  [idl)  and  q>  =  ~(2n/X)  f  A refrdl,  respectively  [11-14].  Using 
the  exit  X-ray  attenuation  and  phase  shift  from  a  given  projection, 
we  then  simulated  the  free-space  Fresnel  diffractions  of  the  exit  X- 
rays  propagating  from  the  sample  to  the  detector  with  a  setting  of 
R2  =  1  m  and  M  =  2  [16].  The  Fresnel  diffraction  intensities  in  all  the 
angular  projections  were  log-transformed  to  calculate  the  ray- 
sums  needed  for  tomographic  image  reconstruction.  We  applied 
the  standard  FDK  algorithm  for  image  reconstruction  [15].  Fig.  1 
shows  a  reconstructed  slice  of  apparent  LAC  values  in  gray  scale. 
This  slice  shows  the  edge-enhancing  dark-bright  fringes  caused 
by  phase  contrast  at  the  interfaces  between  the  spheres.  Fig.  2(a) 
shows  the  profile  of  the  apparent  LAC  values  along  the  vertical  cen¬ 
tral  line  of  Fig.  1.  The  ordinate  denotes  apparent  LAC  value  in  (1/ 
cm)  along  the  vertical  line,  while  the  abscissa  denotes  pixel  loca¬ 
tion  along  the  line.  For  the  sake  of  comparison  Fig.  2(b)  and 
Fig.  2(c)  show  the  profiles  of  the  apparent  LAC  values  along  the  ver¬ 
tical  line,  which  were  computed  by  using  Eq.  (1 )  and  Eq.  (2)  respec¬ 
tively.  Note  that  the  ordinate  scale  in  Fig.  2(b)  is  three  times  larger 
than  that  of  Fig.  2(a)  and  (c).  To  facilitate  the  comparisons  of  the 
three  profiles,  we  measured  the  fringe’s  peak-to-peak  differences 
(the  changes  between  peak  and  trough)  at  eight  fringe  sets  for  each 
of  the  three  profiles  shown  in  Fig.  2(a)-(c).  Table  1  lists  the  mea¬ 
sured  fringe  peak-to-peak  differences  of  apparent  LACs  in 
Fig.  2(a),  which  is  obtained  from  Fig.  1,  the  phantom’s  phase  con¬ 
trast  tomogram  simulated  by  using  the  Fresnel  diffraction  and 
FDK-reconstruction.  Compared  to  these  eight  peak-to-peak  differ¬ 
ences  in  Fig.  2(a),  the  table  lists  as  well  the  errors  of  the  calculated 
peak-peak  differences  based  on  Eq.  (1)  (Fig.  2(b))  and  Eq.  (2) 
(Fig.  2(c))  respectively.  From  these  comparisons,  it  is  clear  that 
the  apparent  LAC  values  calculated  from  Eq.  (1 )  are  grossly  inaccu¬ 
rate  with  the  errors  amounting  to  several  hundred  percent,  while 
the  apparent  LAC  values  calculated  by  means  of  Eq.  (2)  achieve  rea¬ 
sonable  quantitative  accuracies  as  is  shown  in  Table  1.  In  addition, 
the  profile  Fig.  2(c)  shows  two  small  bumps  at  the  center  of  the 
profile.  We  found  that  these  variations  of  apparent  LACs  are  arti¬ 
facts  caused  by  omitting  the  small  nonlinear  terms  in  phase  shifts 
in  the  derivation  of  Eq.  (2).  The  omission  of  these  nonlinear  terms 
mainly  affects  the  accuracies  in  computing  the  side-lobes  of  phase 
contrast  fringes,  thereby  causes  these  two  small  bumps,  which  are 
associated  with  the  distorted  fringe  side-lobes  near  the  right 
boundary  of  the  left  off-center  sphere  in  the  phantom  (Fig.  1). 

In  summary,  we  present  a  general  formula  (Eq.  (2))  to  quantita¬ 
tively  relate  the  apparent  LAC  values  in  cone-beam  phase  contrast 


tomography  to  sample’s  LAC  and  refractive  index  values.  This 
formula  overcomes  the  gross  inaccuracy  of  the  existing  formula 
in  the  literature  in  analyzing  high-resolution  phase  contrast 
tomography.  We  believe  that  this  general  formula  will  be  useful 
for  designing  phase  contrast  tomography  experiments,  and  for  cor¬ 
rectly  interpreting  and  quantifying  the  apparent  LAC  values  in 
cone-beam  X-ray  phase  contrast  tomography.  For  example,  Eq. 
(2)  can  be  used  for  computing  the  edge-enhancement  profiles  for 
analyzing  presampling  modulation  transfer  functions  in  synchro¬ 
tron  radiation  microtomography  [10]. 
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Purpose:  The  robustness  of  the  phase  retrieval  methods  is  of  critical  importance  for  limiting  and 
reducing  radiation  doses  involved  in  x-ray  phase  contrast  imaging.  This  work  is  to  compare  the 
robustness  of  two  phase  retrieval  methods  by  analyzing  the  phase  maps  retrieved  from  the  experi¬ 
mental  images  of  a  phantom. 

Methods:  Two  phase  retrieval  methods  were  compared.  One  method  is  based  on  the  transport  of  in¬ 
tensity  equation  (TIE)  for  phase  contrast  projections,  and  the  TIE-based  method  is  the  most  com¬ 
monly  used  method  for  phase  retrieval  in  the  literature.  The  other  is  the  recently  developed 
attenuation-partition  based  (AP-based)  phase  retrieval  method.  The  authors  applied  these  two  meth¬ 
ods  to  experimental  projection  images  of  an  air-bubble  wrap  phantom  for  retrieving  the  phase  map 
of  the  bubble  wrap.  The  retrieved  phase  maps  obtained  by  using  the  two  methods  are  compared. 

Results:  In  the  wrap’s  phase  map  retrieved  by  using  the  TIE-based  method,  no  bubble  is  recogniz¬ 
able,  hence,  this  method  failed  completely  for  phase  retrieval  from  these  bubble  wrap  images.  Even 
with  the  help  of  the  Tikhonov  regularization,  the  bubbles  are  still  hardly  visible  and  buried  in  the 
cluttered  background  in  the  retrieved  phase  map.  The  retrieved  phase  values  with  this  method  are 
grossly  erroneous.  In  contrast,  in  the  wrap's  phase  map  retrieved  by  using  the  AP-based  method, 
the  bubbles  are  clearly  recovered.  The  retrieved  phase  values  with  the  AP-based  method  are  reason¬ 
ably  close  to  the  estimate  based  on  the  thickness-based  measurement.  The  authors  traced  these  stark 
performance  differences  of  the  two  methods  to  their  different  techniques  employed  to  deal  with  the 
singularity  problem  involved  in  the  phase  retrievals. 

Conclusions:  This  comparison  shows  that  the  conventional  TIE-based  phase  retrieval  method, 
regardless  if  Tikhonov  regularization  is  used  or  not,  is  unstable  against  the  noise  in  the  wrap’s 
projection  images,  while  the  AP-based  phase  retrieval  method  is  shown  in  these  experiments 
to  be  superior  to  the  TIE-based  method  for  the  robustness  in  performing  the  phase  retrieval. 

©2011  American  Association  of  Physicists  in  Medicine.  [DOI:  10.1118/1.3618731] 
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I.  INTRODUCTION 

In  recent  years,  the  in-line  phase-contrast  x-ray  imaging  has 
attracted  intensive  research  efforts.  The  in-line  phase-con¬ 
trast  x-ray  imaging  is  a  technique  using  the  free  space  dif¬ 
fraction  of  phase-shifted  x-rays  to  form  the  interference 
fringes  at  tissues’  boundaries  and  interfaces  in  images.  1-8  In 
fact,  for  over  100  years,  the  tissue  attenuation  differences 
have  been  the  sole  contrast  mechanism  for  medical  x-ray 
imaging.  However,  when  x-rays  traverse  the  body  parts,  as  a 
wave  x-rays  undergo  phase  shifts  as  well.  The  amount  of  the 
phase  shift  along  an  exit  ray  is  determined  by  the  tissue  re¬ 
fractive  indices  along  the  ray.  Note  that  x-ray  refractive 
index  n  is  a  complex  number  and  equal  to  n  =  1  —  8  +  if, 
where  8  is  the  refractive  index  decrement  and  responsible 
for  x-ray  phase  shift,  while  (1  is  the  imaginary  part  of  the  re¬ 
fractive  index  and  responsible  for  x-ray  attenuation.  The 
amount  of  x-ray  phase  shift  along  an  exit  ray  is  given  by 
cp  =  —  (2n/X)  J  Sds,  where  1  is  the  x-ray  wavelength,  and 


the  integral  is  over  the  ray  path. 1  2  In  other  words,  the  phase 
shift  is  equal  to  the  projection  of  refractive  index  decrements 
scaled  by  a  factor  (2n/X).  On  the  other  hand,  x-ray  attenua¬ 
tion  along  an  exit  ray  is  determined  by  the  projection  of  the 
tissue  linear  attenuation  coefficients  along  the  ray,  which  is 
equal  to  (An/ X)  J  fids.1'2  We  have  estimated  8  and  f  values 
for  the  biological  tissues  and  found  that  the  tissue  8  values 
(10~6-10-8)  are  about  1000  times  greater  than  their  j>  values 
(10~9-10-11)  for  10-150  keV  x-rays.8  Hence,  the  differen¬ 
ces  in  x-ray  phase  shifts  between  different  tissues  are  about 
1000  times  greater  than  their  differences  in  the  projected  lin¬ 
ear  attenuation  coefficients.  Therefore,  the  phase-contrast 
imaging  techniques  may  greatly  increase  the  lesion-detection 
sensitivity  for  x-ray  imaging.  The  settings  for  the  inline 
phase-sensitive  imaging  are  similar  to  that  of  conventional 
x-ray  imaging,  provided  a  source  with  very  small  focal  spot 
and  a  sufficiently  large  object-detector  distance  are 
required.5  7'8  In  the  inline  imaging  setting,  x-rays  undergo 
phase  shifts  as  traversing  the  imaged  object,  and  then  diffract 


5073  Med.  Phys.  38  (9),  September  201 1 


0094-2405/201 1/38(9)/5073/8/$$30.00 


©  201 1  Am.  Assoc.  Phys.  Med.  5073 


5074 


Yan,  Wu,  and  Liu:  Robustness  of  phase  retrieval  methods 


5074 


freely  over  a  sufficiently  large  distance  before  reaching  the 
detector.  In  this  way,  the  tissues’  phase  contrast  manifests  as 
the  dark-bright  diffraction  fringes  at  tissues’  boundaries  and 
interfaces  in  the  measured  images.5,7’8  Hence,  the  inline 
phase  contrast  imaging  has  good  potential  of  greatly  enhanc¬ 
ing  the  detection  sensitivity  and  reducing  radiation  doses 
involved  in  the  imaging. 

However,  as  the  interfaces  and  boundaries  of  the  different 
tissue  compartments  are  greatly  accentuated  in  a  phase-con¬ 
trast  image,  the  bulk  phase  contrast  in  a  given  tissue  com¬ 
partment,  where  the  phase  shifts  may  vary  slowly,  could  get 
lost.  This  is  because  the  phase  contrast  is  proportional  to  the 
Laplacian  and  gradient  differentials  of  the  phase  shifts,  as  is 
shown  in  Sec.  II  below.  Moreover,  the  information  about  x- 
ray  phase  shifts  are  valuable  for  tissue  characterizations, 
since  x-ray  phase  shift  along  a  ray  is  proportional  to  the  pro¬ 
jected  tissue  electron  density,  that  is,  cp  —  —Xre  J  peds, 
where  re  is  the  classic  electron  radius  and  pe  denotes  the  tis¬ 
sue  electron  density  and  the  integral  is  over  the  ray  path.1,2 
In  order  to  fully  exhibit  tissue  phase  contrast  and  reconstruct 
tissue  projected  electron  densities  for  quantitative  tissue 
characterizations,  one  needs  to  extract  the  tissue  phase  shifts 
from  the  mixed  contrast  exhibited  in  a  phase-sensitive  pro¬ 
jection.  The  procedure  of  retrieving  the  phase-shift  map  of 
an  object  from  its  phase-sensitive  projections  is  called  the 
phase  retrieval. 10-14  A  retrieved  phase  map  is  able  to  provide 
a  quantitative  map  of  the  object’s  projected  electron  den¬ 
sities,  which  could  be  used  for  quantitative  tissue  character¬ 
izations.12-14  Moreover,  performing  phase  retrieval  is 
necessary  for  reconstructing  volumetric  3-D  maps  of  tissue 
attenuation  coefficients  and  refractive  indices,  respec¬ 
tively,14,15  and  for  eliminating  the  phase-contrast  caused 
artifacts  in  the  volumetric  3-D  images.16 

Examining  the  robustness  of  phase  retrieval  methods 
against  the  projection  noise,  such  as  the  x-ray  photon  quan¬ 
tum  noise  and  others  involved  in  the  projection  acquisitions, 
is  an  important  task  in  phase  contrast  imaging.  If  a  phase  re¬ 
trieval  method  is  not  robust  against  the  noise,  the  method 
could  become  unstable  and  fail  in  phase  retrievals  if  the 
acquired  projection  images  had  be  contaminated  with  sub¬ 
stantial  noise.  On  the  other  hand,  suppressing  x-ray  quantum 
noise  may  require  using  high  radiation  doses  in  the  image 
acquisitions.  For  clinical  applications,  it  is  imperative  to 
limit  and  reduce  radiation  doses  involved,  hence,  it  is  critical 
to  develop  robust  phase  retrieval  methods  for  future  clinical 
applications  of  phase  contrast  imaging.  The  purpose  of  this 
work  is  to  compare  the  robustness  of  two  phase  retrieval 
methods  by  means  of  analyzing  the  phase  maps  retrieved 
from  the  experimental  projection  images  of  an  air-bubble 
wrap  phantom.  We  will  first  introduce  the  phase  retrieval 
method  that  is  based  on  the  transport  of  intensity  equation 
(TIE),  which  describes  how  the  phase  contrast  is  encoded  in 
the  projection  images.10,11  This  TIE-based  method  is  the 
most  commonly  used  method  for  phase  retrieval  in  the  litera¬ 
ture.  We  will  then  introduce  a  recently  developed  phase  re¬ 
trieval  method,  namely  the  attenuation-partition  based  (AP- 
based)  method.13  In  order  to  examine  the  robustness  of  these 
two  methods  against  noises,  we  applied  these  two  phase  re¬ 


trieval  methods  to  experimental  phantom  images,  namely  a 
radiograph  and  a  phase  contrast  image  of  the  air-bubble 
wrap  phantom,  for  the  phase  retrievals.  We  will  analyze  the 
retrieved  phase  maps  obtained  by  using  the  two  methods  and 
compare  their  performances  in  the  robustness  against  noises 
in  the  phantom  images.  In  Sec.  IV,  we  explain  that  these  per¬ 
formance  differences  of  the  two  phase  retrieval  methods  are 
rooted  in  their  different  techniques  employed  to  deal  with 
the  singularity  problem  involved  in  the  phase  retrievals.  In 
addition,  other  factors  that  may  affect  the  phase  retrieval  per¬ 
formance  will  be  briefly  discussed  as  well. 

II.  PHASE  RETRIEVAL  METHODS 

In  phase  contrast  images,  the  attenuation  contrast  and 
phase  contrast  are  mixed  together.  In  order  to  retrieve  phase 
maps  from  phase  contrast  projection  images,  one  should 
understand  how  the  phase  contrast  is  encoded  in  the  projec¬ 
tion  images.  This  understanding  can  be  gained  from  the  x- 
ray  propagation  equation  such  as  the  TIE.10,11  This  equation 
can  be  derived  either  from  the  x-ray  Fresnel-diffraction 
equation  or  the  Wigner  distribution  based  phase-space  for¬ 
malism.4-11  If  we  denote  the  x-ray  transmission  image  of  an 
object,  or  equivalently  its  x-ray  attenuation  map  by  A2(r) 
and  the  x-ray  phase-shift  map  of  the  object  by  cp(r),  then  the 
detected  x-ray  intensities  are  given  by8,1 1,17 

K?d)  =  W  '  (^(^/M)V„(rD/M))  J 

(1) 

where  I  is  the  average  x-ray  wavelength  of  the  polychro¬ 
matic  x-rays.  In  Eq.  (1),  /,„  is  the  entrance  x-ray  intensity, 
R  i  is  the  source-object  distance  set  in  the  projection,  R2  is 
the  object-detector  distance,  M  =  (R\  +  /?2)//?i,  the  magni¬ 
fication  factor  in  the  projection,  and  ip  denotes  the  position 
vector  in  the  detector  plane.  In  Eq.  (1),  the  operator  V 
denotes  the  two-dimensional  transverse  gradient  differential 
operator.  For  the  purpose  of  this  work,  we  assume  that  the 
x-ray  source’s  focus  spot  is  ideally  point-like  and  the  detec¬ 
tor  employed  is  an  ideal  detector.  Obviously,  the  detected 
x-ray  intensity  /(ip)  is  determined  not  only  by  attenuation 
map  A2(r)  in  the  projection,  but  also  by  the  encoded  phase 
contrast,  that  is,  by  the  transverse  Laplacian  and  gradient 
differentials  of  the  phase-shift  map  cp(r)  in  the  projection. 
Since  the  phase  contrast  and  attenuation  contrast  are  mixed 
together  in  a  phase  contrast  image,  as  is  shown  by  Eq.  (1), 
hence  measuring  a  single  phase-contrast  image  is  not 
enough  for  retrieving  the  phase-shift  map  cp(r)  in  the  pro¬ 
jection.  Therefore,  in  general,  at  least  two  projection  images 
are  needed  for  retrieving  the  phase-shift  map  cp(r)  of  the 
object. 

The  most  commonly  used  phase  retrieval  method  in  the  lit¬ 
erature  is  the  TIE-based  method.  In  this  method,  two  acquired 
projection  images  of  an  object  are  used  for  retrieving  the 
phase-shift  map  of  the  object:  one  is  a  contact  radiograph 
A2(r)  of  the  object,  the  other  is  a  phase  contrast  projection 
image  /(ip)  of  the  object.  One  then  is  able  to  retrieve  the 
phase-shift  map  cp(f)  by  solving  Eq.  (1)  for  cp(r)  as  follows12: 
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cp(r)  =  -  (2 7tM//)J?2)V~2 


x 


SVIV- 


(2) 


In  above  equation,  the  operator  V  2  denotes  the  inverse  of 
the  Laplacian  differential  operator  V2  =  (d1  /dx1  +  d2/dy2). 
The  inverse  Laplacian  operator  V  2  is  a  pseudodifferential 
operator.  The  action  on  a  function  g(r )  of  a  pseudodifferential 
operator  such  as  V  2  is  defined  as 


v-2g(f) 


exp(/2jt(5  -  r)  •/)  •  (— 1/^4tt2/*2^ 
x  g(s)d2s-  d2f. 


(3) 


In  Eq.  (3),  s  and /  denote  the  integral  variables  in  the  coordi¬ 
nate  space  and  frequency  space,  respectively.  Hence,  one 
can  compute  the  phase  map  cp(r)  by  using  Eq.  (2)  with  the 
help  of  Eq.  (3).  The  TIE-based  method  is  effective  and  fast 
for  phase  retrievals  in  the  cases  with  strong  phase  contrast 
effects  and  low  noise  levels  as  is  shown  in  many  cases  dis¬ 
cussed  in  the  literature.11'1'  However,  the  Achilles  heel  of 
the  TIE-based  method  lies  at  the  inverse  Laplacian  operator 
V  2  involved  in  Eq.  (2).  The  operator  V  2  adds  a  zero-fre¬ 
quency  singularity  to  Eq.  (2),  as  is  suggested  by  Eq.  (3).  As 
we  will  see  below,  the  singularity  may  amplify  the  noise  ran¬ 
domly  embedded  in  the  acquired  images  and  result  in  insta¬ 
bility  in  phase  retrievals. 

In  order  to  get  rid  of  the  singular  pseudodifferential  oper¬ 
ator  V~2  involved  in  the  TIE-based  phase  retrievals,  we 
recently  developed  a  novel  phase  retrieval  method:  the  AP 
based  method  or  the  AP-based  method  for  short.1"10  Our 
idea  in  this  method  is  to  utilize  the  correlation  between  the 
x-ray  phase  shift  and  its  attenuation  to  eliminate  any  singu¬ 
larity  involved  in  the  phase  retrievals.  As  is  well  known,  in 
the  diagnostic  x-ray  imaging,  x-ray  attenuation  A2  arises 
from  three  x-ray-matter  interactions:  the  photoelectric 
absorption,  the  coherent  scattering,  and  the  incoherent  scat¬ 
tering.  Correspondingly,  we  can  partition  the  x-ray  attenua¬ 
tion  into  two  parts:  the  Compton  scattering-caused 
attenuation  A2KN  and  the  attenuation  caused  by  photoelectric 
absorption  and  coherent  scattering,  which  is  denoted  by 
Ape  coh(r).  That  is,  we  can  partition  the  overall  x-ray  attenua¬ 
tion  A2  into  a  product  of  two  parts 


A20(?)=A2K N(r)-A2pecoh(f). 


(4) 


The  reason  of  factoring  out  A2KN  is  that  both  of  A2KN  and  the 
phase-shift  cp  are  determined  by  the  tissue’s  electron  density 


Am(r)  =  exp 


<p{r)  =  -foepe,p{r)  (5) 


where  re  =  2.818  x  lCL15m  is  the  classical  electron  radius, 
pep ,  the  projected  electron  density  along  the  ray  path,  that  is, 
pep(r)  =  \my  pe{r,z)dz,  and  Gkn  denotes  the  average  Comp¬ 
ton  scattering  cross-section  for  polychromatic  x-rays.  Note 
that  Compton  scattering  cross-section  is  given  by  the  Klein- 
Nishina  total  cross-section19 


Okn(E)  =2nr2i]-^-  '■ 

l  T  L 


2(1  +  //) 


1  +  2>i 


-log(l  +2  rj) 
1 1 


+  — log(l  +  2f/) 

2/7 


1  +  3i/  ) 

(1  +  2)/)2  J  ’ 


(6) 


where  ij  =  E  /5 1 1  keV  and  it  slowly  varies  with  x-ray 
energy  E  for  E<C  51 1  keV.  We  observed  that  the  extent  of 
the  correlation  between  phase  and  attenuation  depends  on 
the  x-ray  photon  energies  and  the  tissues  elemental  composi¬ 
tions.  For  example,  for  light  elements  with  atomic  numbers 
Z  <  10  and  x-rays  of  60  keV  or  higher,  the  x-ray-matter 
interactions  are  dominated  by  the  Compton  scattering, 
hence,  both  the  tissue  attenuation  and  phase  shift  are  all 
determined  by  tissues’  electron  density  distributions.  We  call 
this  relationship  between  phase  shift  and  attenuation  the 
phase-attenuation  duality,  and  we  define  the  so-called  duality 
transform  as20 

t>(/(rD))  =  (1  -  (I2 R2re/2nMdKN)  V2)  • 

(7) 

For  the  cases  where  the  x-ray-matter  interactions  are  domi¬ 
nated  by  the  Compton  scattering,  we  proved  that  the  attenua¬ 
tion  A2kn  can  be  found  from  a  single  phase  contrast  image 
1(?d)  by  means  of  the  duality  transform  such  that  A2KN{r) 
=  ®(/(rb)).20  With  the  help  of  the  phase-attenuation  duality, 
the  phase-shift  map  of  the  object  can  be  simply  retrieved  as 

cp(r)  =  —In  A2KN(f)  =  ^-lnX>(/(fb)).  (8) 


We  call  this  formula  as  the  phase-attenuation  duality-based 
phase  retrieval  formula.20  One  significant  advantage  of  this 
duality-based  phase  retrieval  formula  is  that  the  pseudodif¬ 
ferential  operator  (l  —  (2?R2re/2nMdKN)\72)  involved  in 
Eqs.  (7)  and  (8)  is  free  of  any  singularity.  Therefore,  the 
phase  retrieval  method  based  on  the  phase-attenuation  dual¬ 
ity  is  stable  and  robust  against  the  noise  in  images.20 

Unfortunately,  the  phase-attenuation  duality  does  not  hold 
in  some  other  cases  such  as  imaging  with  low  energy  x-rays  or 
imaging  bones  and  calcifications  that  contain  substantial 
amount  of  heavy  elements.  In  those  cases,  the  retrieved  phase 
map  using  Eq.  (8)  may  provide  merely  coarse  estimates  of  the 
real  phase  maps  of  the  imaged  objects.  In  order  to  overcome 
this  limitation  of  the  duality-based  phase  retrieval  formula  Eq. 
(8),  recently  we  proposed  to  correct  the  errors  iteratively 
through  repeated  comparisons  of  the  computed  estimates  of 
the  x-ray  Fresnel-diffraction  intensities  against  the  acquired 
projection  intensities.  To  implement  this  strategy,  we  devel¬ 
oped  the  attenuation-partition  based  iterative  algorithm,  whose 
flow  chart  is  shown  in  Fig.  I.13  While  interested  readers  can 
find  the  mathematical  proofs  of  the  algorithm  in  Ref.  13,  here, 
we  give  a  brief  explanation  of  this  algorithm  flow  chart  in 
Fig.  1.  In  this  flow  chart.  A2  denotes  the  attenuation  map  of 
the  imaged  object  measured  from  its  contact  radiograph  and  I 
is  the  acquired  phase  contrast  image  of  the  object.  In  the  flow 
chart,  ! kn  denotes  the  simulated  phase  contrast  image  formed 
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Fig.  1 .  Flow  chart  of  the  AP  based  iterative  algorithm. 

by  using  the  estimated  object’s  electron  densities  and  by 
turning  off  the  attenuation  processes  associated  with  x-ray 
photoelectric  absorption  and  coherent  scattering.  In  Fig.  1,  two 
transforms  are  implemented  in  the  iterations,  one  is  the  Fres- 
nel-diffraction  transform  grc>  which  implements  x-ray 
Fresnel-diffraction,21  and  the  other  is  the  duality  transform  Xs 
as  is  defined  in  Eq.  (7).  According  to  this  flow  chart  (Fig.  1)  of 
the  AP-based  phase  retrieval  method,  one  first  applies  the 
duality  transform  T>  to  the  acquired  phase  contrast  image  / 
to  obtain  an  estimate  of  A2KN  and  the  phase  map  cp  of  the 
object.  One  then  compares  A2KN  against  the  measured  A 2 
and  computes  the  error  8 A  =  —  Aq.  Applying  the  Fres¬ 

nel-diffraction  transform  g  t0  the  fictitious  transmission 
function  8A  ■  exp (i<p),  one  computes  the  I kn  -correction 
81  =  |5re(M  •  exp(i<p))|  ,  and  a  new  estimate  of  I ^  is  deter¬ 
mined  as  7/av  =  {V~l  +  V8I)1-  With  this  new  estimate  of  1^, 
one  can  start  a  new  round  of  the  iteration  as  is  indicated  in  the 
flow  chart  Fig.  1.  The  iteration  converges  when  the  IKN  does 
not  change  substantially  with  further  iteration  steps.1 1 

III.  COMPARISON  OF  THE  TWO  PHASE  RETRIEVAL 
METHODS 

In  order  to  compare  the  robustness  of  above  two  phase  re¬ 
trieval  methods,  namely,  the  TIE-based  method  and  the  AP- 
based  method,  we  applied  the  two  methods  for  retrieving  the 
phase  map  of  a  piece  of  air-bubble  wrap.  In  the  bubble  wrap, 
the  air  packets  are  locked  between  two  thin  layers  of  low-den¬ 
sity  polyethylene  films  forming  the  air  bubbles.  Shown  in  Fig. 
2(a)  is  a  contact  radiograph  of  the  wrap.  The  x-ray  source 
used  was  a  micro-focus  x-ray  tube  (L8 121-02,  Hamamatsu) 
with  a  focal  spot  size  of  7  /mi.  The  source  has  a  tungsten  tar¬ 
get  without  added  filtration  and  was  operating  at  40  kVp  and 
0.2  mA  for  a  30  s  exposure.  The  average  x-ray  photon  energy 
was  estimated  at  14.6  keV  for  the  40  kVp  x-rays.  In  this  ac¬ 
quisition,  the  source-detector  distance  (SID)  was  set  to  1.78 
m.  The  imaging  detector  employed  was  an  aSe-based  flat- 
panel  detector  (DirectRay,  Hologic)  with  a  pixel  pitch  of  140 
/mi.  Since  the  thin  polyethylene  films  and  locked  air  bubbles 
in  the  wrap  generate  little  differences  in  x-ray  attenuation, 
hence,  the  wrap’s  radiograph  in  Fig.  2(a)  exhibits  only  little 
image  contrast,  and  the  x-ray  quantum  noise  in  the  acquisition 
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Fig.  2.  (a)  Contact  radiograph  of  the  bubble  wrap  acquired  at  40  kVp  and 
with  a  SID  =  1.75  m.  (b)  Intensity  profile  along  the  marked  line  on  (a). 

masked  the  wrap’s  details  such  that  the  rims  of  many  bubbles 
are  not  visible  on  the  image.  Figure  2(b)  shows  a  profile  of 
image  intensities  along  the  marked  dash  line  on  Fig.  2(a)  and 
the  profile  reveals  how  the  noisy  intensity  variations  mask  the 
low  contrast  rims  of  the  bubbles.  In  contrast.  Fig.  3(a)  is  a 
phase  contrast  image  of  the  bubble  wrap.  This  phase  contrast 
image  was  acquired  with  a  sample-detector  distance 
R2  =  1.15  m  and  otherwise  the  same  technique  settings  as 
that  used  for  the  radiograph,  that  is,  with  the  same  focal  spot 
size,  same  tube  voltage,  same  tube  current  and  exposure  time, 
same  SID,  and  the  same  detector.  However,  in  this  projection, 
the  phase-shifted  x-rays  were  allowed  to  freely  diffract  over 
1.1 5-m  long  distance  on  their  way  to  the  detector  and  to  form 
the  interference  fringes  at  the  bubble  boundaries  on  the  image. 
As  is  shown  in  Fig.  3(a),  the  phase  contrast  depicts  not  only 
the  bubble  rims,  but  also  the  dents  inside  the  bubble  domes. 
The  enhancement  is  resulted  from  the  rapid  changes  of  the 
projected  thicknesses  of  the  polyethylene  films  at  the  bubble 
rims  and  the  dents  inside  bubble  domes.  These  rapid  changes 
in  the  projected  thicknesses  of  polyethylene  generate  as  well 
rapid  changes  in  x-ray  phase  shifts,  since  the  x-ray  phase 
shifts  are  given  by  <p(f)  =  — XrepcTp{f ),  where  pe  denotes 
polyethylene  film’s  electron  density  and  Tp(r)  denotes  pro¬ 
jected  thickness  of  polyethylene  along  the  ray.  According  to 
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(a) 


(b) 


Fig.  3.  (a)  Phase  contrast  image  of  the  bubble  wrap,  which  was  acquired  at 
40  kVp  and  with  a  SID  =1.75  m  and  a  magnification  factor  M  =  2.8. 
(b)  Intensity  profile  along  the  marked  line  on  Fig.  3(a). 


Eq.  (1),  the  Laplacian  and  gradient  differentials  of  the  rapid 
phase  shifts  generated  large  x-ray  intensity  variations,  which 
manifest  as  the  bright-dark  fringes  flank  along  the  bubble  rims 
in  Fig.  3(a).  Figure  3(b)  shows  a  profile  of  image  intensities 
along  the  marked  dash  line  on  Fig.  3(a),  as  the  line  traces  the 
same  positions  as  indicated  by  the  dash  line  on  Fig.  2(a).  The 
downward  and  upward  overshooting  of  the  intensity  values  in 
the  profile  represent  the  bright  and  dark  fringes  at  the  bubble 
boundaries  in  phase  contrast  image  Fig.  3(a).  These  up-down 
bipolar  overshootings  of  intensities  are  much  larger  than  the 
background  noise,  hence,  the  noise  gets  masked  in  Fig.  3(a). 
While  the  edge-enhancement  generated  by  phase-contrast  is 
generally  useful  for  imaging  the  wrap,  however,  such  edge- 
enhancements  may  lead  interpretation  errors  in  the  characteri¬ 
zation  of  the  wrap’s  structure  and  composition.  For  example, 
in  Fig.  3(a),  the  dark  fringes  flanked  the  enhanced  bubble  rims 
may  be  confused  with  possible  structural  breaks  in  the  wrap 
sample,  although  we  knew  in  advance  that  this  wrap  sample  is 
free  of  any  damage.  In  order  to  fully  exhibit  the  wrap’s  phase 
contrast  and  correctly  characterize  the  bubble  wrap,  it  is  nec¬ 
essary  to  perform  the  phase  retrieval.  In  order  to  retrieve  the 
wrap’s  phase-shift  map  from  the  wrap’s  radiograph  [Fig.  2(a)] 
and  its  phase  contrast  image  [Fig.  3(a)],  we  first  employed 
the  TIE -based  phase  retrieval  formula  Eq.  (2).  The  retrieved 


phase  map  with  the  TIE-method  is  shown  Fig.  4(a).  Appa¬ 
rently,  the  phase  retrieval  failed  completely,  since  no  bubble 
is  recognizable  in  Fig.  4(a).  Although  the  failure  can  be 


Fig.  4.  Retrieved  phase  maps  of  the  bubble  wrap,  (a)  With  the  TIE-based 
method,  (b)  With  the  TIE-based  method  and  Tikhonov  regularization,  (c) 
Profile  of  the  retrieved  phase  values  along  the  marked  line  on  (b). 
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ascribed  to  the  high-levels  of  the  quantum  noise  presented  in 
the  radiograph  of  Fig.  2(a),  but  it  is  rooted  in  the  intrinsic 
instability  of  the  TIE-method.  As  we  pointed  out  in  Sec.  II, 
the  inverse  Laplacian  operator  V  2  plugs  in  the  zero-fre¬ 
quency  singularity  into  the  TIE-base  phase  retrieval  formula 
Eq.  (2).  The  singularity  amplified  the  noise  randomly  embed¬ 
ded  in  the  projections  and  results  in  the  failure  in  the  phase 
retrievals.  To  mend  the  instability  of  the  TIE-method,  one  can 
try  to  apply  the  regularization  techniques  to  tame  the  singular¬ 
ity-caused  problems.  A  common  regularization  technique  in 
the  literature  is  Tikhonov  regularization.22  With  this  regulari¬ 
zation  scheme,  one  replaces  the  inverse  Laplacian  operator 
V  2  by  a  pseudodifferential  operator  V2/[(V2)2  +  a2],  where 
a  is  the  regularization  parameter,  which  is  roughly  propor¬ 
tional  to  the  images  noise-signal  ratios  of  the  images.22  In 
essence,  Tikhonov  regularization  seeks  the  minimum-norm, 
least  squares  solution  of  Eq.  (2).  While  the  Tikhonov  regulari¬ 
zation  of  the  inverse  Laplacian  operator  may  tame  the  insta¬ 
bility,  it  sacrifices  the  phase  retrieval  accuracies.  The 
performance  in  phase  retrieval  of  Tikhonov  regularization  is 
highly  dependent  on  the  amounts  of  noise  presented  in  the 
images.  The  regularization  parameter  a  was  selected  through 
a  trial  by  comparing  the  phase  maps  retrieved  with  a  wide 
range  of  a-values  such  as  a  =  Ah2,  5  Ah2,  10  Ah2,  50  Ah2, 
100  Ah2,  200  Ah2,  300  Ah2,  and  500  Ah2,  where  Ah  is  the  fre¬ 
quency  sampling-step  used  in  the  phase  retrieval.  Figure  4(b) 
is  the  wrap  phase  map  retrieved  using  the  TIE-method  and 
Tikhonov  regularization  with  a  =  200  and  Ah2  =  3.782 
x  1 0  s  /am  2 ,  which  was  the  a  -value  for  the  best  results  as 
determined  from  the  trial.  Apparently,  these  retrieval  results 
are  very  unsatisfactory.  In  Fig.  4(b),  the  bubble  rims  are 
hardly  visible  in  the  extremely  cluttered  background.  Figure 
4(c)  shows  the  profile  of  the  retrieved  phase  values  along  the 
marked  dash  line  on  Fig.  4(b).  As  we  will  explain  below  that 
these  phase  values  are  grossly  erroneous. 

As  a  comparison,  we  have  applied  the  AP-based  method 
to  the  wrap’s  images  in  Figs.  2(a)  and  3(a)  for  the  phase  re¬ 
trieval  as  well.  Following  the  method’s  iteration  flow  chart 
in  Fig.  1,  we  succeeded  in  retrieving  the  wrap’s  phase  map 
of  the  bubble  wrap  with  ten  iteration  steps.  Figure  5(a)  is  the 
retrieved  phase  map  of  the  wrap  with  the  AP-based  method. 

In  a  stark  contrast  to  Figs.  4(a)  and  4(b)  where  the  bubble 
rims  are  hardly  visible,  the  bubble  rims  are  prominently 
depicted  in  Fig.  5(a).  A  profile  of  retrieved  phase  values 
along  the  marked  line  on  Fig.  5(a)  is  shown  in  Fig.  5(b). 
Note  that  all  the  marked  lines  in  Figs.  3(a),  4(b),  and  5(a), 
respectively  correspond  to  the  same  line  defined  on  the  bub¬ 
ble  wrap.  In  this  way,  one  can  easily  compare  the  three  pro¬ 
files  shown  in  Figs.  3(b),  4(c),  and  5(b).  In  order  to  study  the 
quantitative  aspects  of  the  retrieved  phase  maps  of  the  wrap, 
remember  that  the  amount  of  phase  shift  along  a  ray  is  given 
by  cp(f)  =  —ArepeTe(r),  here,  pe  denotes  polyethylene  film’s 
electron  density  and  Tp(r)  denotes  projected  thickness  of 
polyethylene  along  the  ray.  Note  that  x-ray  phase  shifts 
should  be  of  negative  values,  as  x-ray  refractive  indices  of 
tissues  and  materials  are  complex  and  their  real  part  are  less 
than  one.  The  x-rays  traversed  longer  paths  in  polyethylene 
at  the  bubble  rims  and  incurred  larger  phase  shifts,  as 
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Fig.  5.  (a)  Retrieved  phase  maps  of  the  bubble  wrap  with  the  AP-based 
method,  (b)  Profile  of  the  retrieved  phase  values  along  the  marked  line  on  (a). 


compared  to  other  parts  of  the  wrap.  For  example,  the  three 
sharp  negative  peaks  in  the  profile  Fig.  5(b)  reflect  the  large 
x-ray  phase  shifts  incurred  at  the  three  rim  locations  along 
the  marked  line  in  Fig.  5(a).  In  order  to  gauge  the  accuracies 
of  the  retrieved  phase  values  profiled  in  Fig.  5(b),  one  has  to 
know  the  values  of  Tp(r),  the  projected  thicknesses  in  poly¬ 
ethylene  films.  Apparently,  Tp(r)  depends  on  the  local  curva¬ 
tures  of  bubble  domes  and  the  incident  angles  of  the  rays. 
While  it  is  hard  to  measure  Tp(f)  values  in  the  bubbles,  we 
measured  the  thicknesses  of  the  flat  bases  between  bubbles 
using  a  caliber  ruler,  and  we  found  that  the  base  thickness 
was  about  0.055  mm.  Knowing  that  the  molecular  formula 
of  the  low-density  polyethylene  film  is  (C2H4)n  and  its  mass 
density  is  0.925  g/cm3,  we  calculated  its  electron  density  as 
3.18  x  1023/cm\  Assuming  approximately  all  rays  have 
normal  incidence,  we  find  that  the  approximate  phase  shifts 
at  the  flat  bases  between  bubbles  is  about  —4.2  radian.  On 
the  other  hand,  according  to  the  phase  profile  in  Fig.  5(b), 
the  average  phase  shift  at  the  flat  bases  between  bubbles  is 
—  1.05  radian.  Obviously,  these  two  estimates  of  the  phase 
values  though  different,  but  are  reasonably  close.  In  compar¬ 
ison,  the  phase  profile  in  Fig.  4(c),  which  was  obtained  by 
using  the  TIE-based  method  with  Tikhonov  regularization, 
depicts  just  messy  up-down  peaks  buried  in  cluttered  and 
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noisy  background.  In  the  phase  profile.  Fig.  4(c)  only  two 
negative  peaks  at  the  rims  are  barely  recognizable,  and  the 
third  expected  negative  peak  completely  disappears.  Accord¬ 
ing  to  this  profile,  the  phase-shift  values  at  the  bases  fluctu¬ 
ate  over  a  wide  range  from  0  to  +87.5  radian  with  an 
average  of  +27.3  radian,  while  x-ray  phase  shifts  should  be 
negative.  Hence,  the  average  phase  value  of  flat  bases  esti¬ 
mated  based  on  the  profile.  Fig.  4(c)  is  in  a  gross  discrepancy 
with  the  thickness  measurement-based  estimate.  Therefore, 
the  above  comparisons  show  clearly  that  the  AP-based  phase 
retrieval  method  is  shown  to  be  superior  to  the  TIE-based 
method,  regardless  if  the  Tikhonov  regularization  is  used,  in 
performing  the  phase  retrieval  for  the  bubble  wrap. 

Before  ending  this  section,  we  want  to  demonstrate  the  use¬ 
fulness  of  performing  the  phase  retrieval  for  tissue  and  mate¬ 
rial  characterizations  by  comparing  the  bubble  wrap’s  phase 
map  in  Fig.  5(a)  to  its  phase  contrast  image  in  Fig.  3(a).  As  we 
pointed  out  earlier,  the  dark  fringes  which  are  flanking  the 
bubble  rims  in  the  phase  contrast  image  Fig.  3(a),  may  be  con¬ 
fused  for  indicating  possible  structural  breaks  in  the  bubble 
wrap.  In  contrast  to  Fig.  3(a)  riddled  with  the  dark  fringes,  the 
retrieved  phase  map  Fig.  5(a)  does  not  present  any  dark  fringes 
at  all.  This  observation  is  also  verified  by  comparing  the  phase 
profile  Fig.  5(b)  to  the  intensity  profile  Fig.  3(b).  The  up-down 
bipolar  overshootings  appearing  in  Fig.  3(b)  disappear  in  the 
phase  profile  Fig.  5(b),  where  only  sharp  downward  peaks 
present  for  depicting  the  large  phase  shifts  at  the  bubble  rims. 
Therefore,  the  wrap’s  phase  map  Fig.  5(a)  clarifies  the  nature 
of  the  dark  fringes  in  Fig.  3(a)  such  that  these  dark  fringes  do 
not  indicate  any  structural  break  in  the  wrap  sample. 

IV.  DISCUSSION  AND  CONCLUSIONS 

This  study  demonstrates  that  the  robustness  against  noises 
of  a  phase  retrieval  method  is  critical  for  the  qualities  of  the 
retrieved  phase  map.  Specifically,  the  striking  differences 
between  the  retrieved  phase  maps  in  Figs.  4(a)  and  4(b), 
and  Fig.  5(a)  underline  the  importance  of  removing  the 
zero-frequency  singularity  that  is  intrinsic  to  the  TIE-based 
method.  While  the  TIE-based  phase  retrieval  method  is  compu¬ 
tational  effective  for  the  cases  with  strong  phase  contrast  effects 
and  low  noise  levels,11’12  it  completely  failed  in  retrieving  the 
phase  map  of  the  bubble  wrap,  as  is  shown  in  Figs.  4(a)  and 
4(b).  This  is  because  the  noise  problem  is  especially  challeng¬ 
ing  for  the  case  at  hand,  as  the  bubble  wrap  has  very  low 
attenuation  contrast  and  only  moderate  x-ray  exposures  were 
applied  in  bubble  image  acquisitions.  From  mathematical  view¬ 
point,  the  TIE-based  phase  retrieval  formula  Eq.  (2)  is  ill-posed, 
since  it  involves  the  inverse  Laplacian  operator  V  2  that  owns 
a  zero-frequency  singularity.  Beyond  the  mathematics  formal¬ 
ity,  this  zero-frequency  singularity  reflects  the  fact  that  the 
phase  contrast  projection  is  insensitive  to  the  slow  variations  of 
x-ray  phase  shifts.  Seeking  to  recover  the  phase  shifts  by  com¬ 
paring  the  two  projection  images  contaminated  with  noise,  the 
TIE-based  method  has  amplified  the  noise  in  retrieving  the 
slowly  varying  phase  shifts  and  ruined  the  phase  retrieval  com¬ 
pletely.  The  Tikhonov  regularization  mends  the  singularity 
problem  by  replacing  the  inverse  Laplacian  operator  V  2  in 


Eq.  (2)  by  a  pseudodifferential  operator  V2/[(V2)2+oc2]  with  a 
regularization  parameter  a2  >  0.  As  is  shown  in  Fig.  4(b),  this 
regularization  made  some  but  only  little  improvement  com¬ 
pared  to  Fig.  4(a)  as  the  bubbles  start  to  appear  in  the  phase 
map  but  are  still  hardly  visible.  From  the  mathematical  formu¬ 
lation  of  the  Tikhonov  regularization  technique,  it  is  clear  that 
Tikhonov  regularization  just  limits  the  overall  magnitude  of  the 
errors  caused  by  the  noise  in  the  retrieved  phase  map,  but  the 
details  of  the  original  phase  map  could  still  get  lost  in  the  re¬ 
trieval.2"  The  AP-based  method  mends  the  singularity  problem 
by  utilizing  the  correlations  between  x-ray  attenuation  and  x- 
ray  phase  shifts  generated  by  tissues  or  materials,  as  is  shown 
in  Eq.  (8).  In  a  sense,  the  AP-based  method  utilizes  the  object’s 
x-ray  attenuation  to  get  rid  of  the  low-frequency  singularity  and 
the  associated  noise  amplification.  This  physics-motivated  reg¬ 
ularization  scheme  for  the  AP-based  method  is  expected  to  be 
more  effective  than  Tikhonov  regularization  in  phase  retrievals. 
The  greatly  improved  quality  of  the  retrieved  phase  map  shown 
in  Fig.  5(a),  as  compared  to  Figs.  4(a)  and  4(b),  validates  this 
expectation.  Hence,  the  performance  differences  of  these  phase 
retrieval  methods  underscore  the  importance  of  developing  an 
effective  method  to  remove  the  singularity  that  is  intrinsic  to 
the  TIE-based  method.  After  all,  the  high  noise  levels  in  the 
acquired  images  still  take  tolls  on  the  quality  of  the  phase  map 
Fig.  5(a)  retrieved  the  AP-method,  where  the  noise  are  quite 
visible.  While  increasing  radiation  doses  used  in  the  projections 
is  a  possible  solution  for  improving  phase  retrievals,  but  the 
radiation  dose  constraints  are  stringent  in  many  applications.  In 
clinical  applications,  it  is  imperative  to  limit  and  reduce  radia¬ 
tion  doses  involved  in  the  imaging.  Therefore,  more  research  is 
needed  for  developing  ever-improved  phase  retrieval  methods 
for  future  clinical  applications  of  phase  contrast  imaging. 

Other  factors,  which  are  less  critical  to  the  robustness  of 
phase  retrieval  but  may  affect  the  accuracies  of  retrieved 
phase  maps,  include  the  detector  calibrations  such  as  the  flat- 
field  and  gain  corrections,  and  the  ways  to  incorporate  the 
spectral  averaging  effects  of  polychromatic  x-rays.  In  fact, 
the  residual  background  nonuniformity  in  the  acquired 
images  [Figs.  2(a)  and  3(a)]  also  caused  variations  in  the 
background  phase  values  in  Fig.  5  and  reduced  their  accura¬ 
cies.  A  careful  detector  calibration  in  future  experiments 
should  avoid  this  kind  of  problem.  In  addition,  note  that  x- 
ray  wavelength  enters  as  an  important  parameter  in  Eq.  (2) 
for  the  TIE-based  method,  and  in  the  flow  chart  in  Fig.  1  for 
the  AP-based  method.  The  average  x-ray  wavelength  used  in 
the  phase  retrievals  should  represent  the  wavelength’s  linear 
and  nonliner  effects  averaged  over  the  exiting  x-rays  spec¬ 
trum  and  the  detector’s  spectral  response.  For  the  TIE-based 
method,  several  works  discussed  the  ways  of  performing  the 
spectral  averages  over  polychromatic  x-rays.8,23  We  did  not 
use  these  techniques  in  this  work,  because  the  TIE-based 
method  in  any  way  failed  the  wrap’s  phase  retrieval.  In  this 
work,  we  simply  use  the  estimated  average  photon  energy  of 
the  incident  x-ray  to  incorporating  the  spectral  averaging.  In 
the  future,  we  may  develop  more  elaborated  ways  for  incor¬ 
porating  the  spectral  averaging  effects.  Finally,  we  mention 
that  one  could  as  well  employ  a  special  phase  retrieval 
method  for  the  simple  samples  such  as  the  bubble  wrap.  This 
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special  method  works  specifically  for  the  single-material  ho¬ 
mogeneous  samples,  as  long  as  the  linear  attenuation  coeffi¬ 
cient  and  refractive  index  of  this  material  are  provided.  In 
this  special  method,  the  phase  map  of  a  single-material  sam¬ 
ple  can  be  retrieved  from  just  a  single  phase  contrast  image 
of  the  sample.24  The  air-bubble  wrap  can  be  approximately 
treated  as  a  single-material  sample  as  the  air  in  bubbles  con¬ 
tributes  negligibly  to  x-ray  attenuation  and  phase  shifts.  In 
this  work,  we  do  not  compare  this  special  method  for  single¬ 
material  samples  to  the  TIE-based  and  AP-based  methods, 
since  both  the  TIE-based  and  AP-based  methods  are  the  gen¬ 
eral  phase  retrieval  methods  applicable  for  any  general 
objects  of  inhomogeneous  materials. 

In  summary,  in  this  work,  we  compared  the  robustness  of 
two  phase  retrieval  methods,  the  TIE-based  method  and  the 
AP-based  method,  by  analyzing  the  retrieved  phase  maps 
from  the  experimental  projection  images  of  an  air-bubble 
wrap.  We  showed  that  the  TIE-based  method,  regardless  if 
the  Tikhonov  regularization  is  used,  failed  in  retrieving  the 
wrap’s  phase  map.  In  contrast,  in  the  wrap  phase  map 
retrieved  by  using  the  AP-based  method  bubbles  are  clearly 
recovered.  The  retrieved  phase  values  with  this  method  are 
reasonably  close  to  the  estimate  based  on  the  thickness- 
based  measurement.  The  stark  performance  differences  of 
the  two  methods  are  rooted  in  their  different  techniques 
employed  to  deal  with  the  singularity  problem.  This  compar¬ 
ison  shows  that  the  conventional  TIE-based  phase  retrieval 
method,  regardless  of  using  Tikhonov  regularization  or  not, 
is  unstable  against  the  noise  in  the  wrap’s  projection  images, 
while  the  AP-based  phase  retrieval  method  is  shown  in  these 
experiments  to  be  superior  to  the  TIE-based  method  for  the 
robustness  in  performing  the  phase  retrieval. 
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Abstract.  The  goal  of  this  preliminary  study  was  to  investigate  the  effects  of  x-ray  beam  hardening  on  the  detective  quantum 
efficiency  (DQE)  and  the  radiation  dose  of  an  inline  x-ray  imaging  system.  The  ability  to  decrease  the  risk  of  harmful  radiation 
to  the  patient  without  compromising  the  detection  capability  would  more  effectively  balance  the  tradeoff  between  image  quality 
and  radiation  dose,  and  therefore  benefit  the  fields  of  diagnostic  x-ray  imaging,  especially  mammography.  The  DQE  and  the 
average  glandular  dose  were  both  calculated  under  the  same  experimental  conditions  for  a  range  of  beam  hardening  levels, 
corresponding  to  no  added  beam  hardening  and  two  thicknesses  each  of  Rhodium  (Rh)  and  Molybdenum  (Mo)  filters.  The 
dose  calculation  results  demonstrate  a  reduction  of  15%  to  24%  for  the  range  of  beam  hardening  levels.  The  comparison  of 
all  quantities  comprising  the  DQE  exhibit  very  close  correlation  between  the  results  obtained  without  added  beam  hardening 
to  the  results  corresponding  to  the  range  of  beam  hardening  levels.  For  the  specific  experimental  conditions  utilized  in  this 
preliminary  study,  the  results  are  an  indication  that  the  use  of  beam  hardening  holds  the  potential  to  reduce  the  radiation  dose 
without  decreasing  the  performance  of  the  system.  Future  studies  will  seek  to  apply  this  method  in  a  clinical  environment  and 
perform  a  comprehensive  image  quality  evaluation,  in  an  effort  to  further  evaluate  the  potential  of  beam  hardening  to  balance 
the  tradeoff  between  dose  and  image  quality. 

Keywords:  Diagnostic  x-ray  imaging,  beam  hardening,  radiation  dose,  detective  quantum  efficiency  (DQE),  noise  equivalent 
quanta  (NEQ),  modulation  transfer  function  (MTF),  noise  power  spectrum  (NPS) 


1.  Introduction 

Throughout  the  history  of  mammography,  the  technology  and  methods  have  evolved  with  a  consistent 
goal  of  improving  the  tradeoff  between  the  risk  of  harmful  radiation  dose  to  the  patient  and  the  benefit 
of  disease  detection  provided  by  improved  image  quality  [1-4].  However,  the  physical  formulation  of 
the  images,  which  relies  solely  on  attenuation  contrast,  has  remained  the  same.  The  similar  composition 
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of  benign  and  malignant  breast  tissue  produces  very  low  attenuation  contrast  [5-7],  which  presents  a 
considerable  challenge  for  cancer  detection  in  mammography.  Unfortunately,  current  techniques  utilized 
in  mammography  to  improve  the  image  quality  also  deliver  an  increased  radiation  dose.  As  a  result,  the 
challenge  remains  to  more  effectively  optimize  the  tradeoff  between  image  quality  and  radiation  dose. 


2.  Background 

Breast  cancer  typically  arises  in  the  glandular  tissue  [5,8];  therefore,  the  average  glandular  dose  has 
been  established  as  the  standard  measurement  of  radiation  dose  in  mammography,  and  guidelines  have 
been  created  by  numerous  national  and  international  councils  for  its  calculation  and  supervision  in  clinical 
environments  [8-12].  The  formula  for  calculation  of  the  average  glandular  dose  Dg  is  as  follows  [5,13, 
14]: 

Dg  =  DgN  ■  XEse,  (1) 

where  T)gyr  is  the  normalized  average  glandular  dose  coefficient  and  Xese  is  the  object  entrance 
exposure.  I)gjy  is  determined  by  experimental  and  computer  simulation  methods  based  on  the  following 
factors:  radiation  quality  (x-ray  energy  or  half  value  layer),  x-ray  tube  target  material,  filter  material, 
breast  thickness  and  breast  tissue  composition  [5].  Due  to  the  complicated  nature  of  the  calculations, 
as  well  as  the  small  range  of  values  in  mammography  for  each  of  the  T)gy-  calculation  parameters, 
numerous  studies  have  published  tables  of  Dgy-  values  that  are  widely  used  in  clinical  and  research 
environments  [13-16].  However,  the  values  must  typically  be  calculated  directly  in  investigational 
studies,  due  to  the  use  of  parameters  outside  the  standard  values.  This  study  will  therefore  apply  a  Monte 
Carlo  process  described  previously  [14,17,18]. 

A  comprehensive  quantitative  determination  of  the  performance  of  an  x-ray  imaging  system  is  the 
detective  quantum  efficiency  (DQE)  [5,19,20],  as  it  represents  the  transfer  function  of  the  SNR  of  the 
system,  which  is  demonstrated  by  the  following  equation  [5,19,21-24]: 

SNRlUT  S(0)2  ■  MTF(f)?  NEQ(f) 

DQE(,)  =  -SNRff  =  NPSU)  ■  0  ® 

where  the  modulation  transfer  function  (MTF)  is  a  measurement  of  the  frequency  response  of  the 
system  [5,19,23,25,26],  the  noise  power  spectrum  (NPS)  is  a  determination  of  the  noise  processed  by 
the  system  [5,19,23,27-29],  and  5(0)  is  the  mean  intensity  of  a  large  area  signal,  assuming  a  uniform 
x-ray  beam  absorbed  by  the  detector.  The  noise  equivalent  quanta  (NEQ)  is  a  combination  of  the 
three  quantities  and  describes  the  number  of  quanta  recorded  at  each  spatial  frequency,  which  provides 
information  regarding  the  maximum  attainable  value  of  the  output  SNR  [22,30-33].  The  photon  fluence 
( q )  is  the  number  of  incident  photons  per  square  millimeter,  which  is  equal  to  the  square  of  the  input  SNR 
for  x-ray  photons  following  a  Poisson  distribution  [5].  In  this  preliminary  study,  comparisons  were  made 
of  the  DQE  as  well  as  the  individual  quantities,  in  an  effort  to  provide  a  comprehensive  investigation  of 
the  effects  of  beam  hardening  on  the  imaging  performance. 

Beam  hardening  is  the  removal  of  low  energy  photons  from  the  x-ray  beam  before  reaching  the  object, 
through  the  use  of  filtration  or  other  objects  placed  in  the  path  of  the  x-ray  beam  [5].  Low  energy  photons 
are  absorbed  more  readily  due  to  lower  penetration  ability  [5] ;  therefore,  removing  them  prior  to  x-ray 
exposure  has  the  potential  to  reduce  the  patient  dose  without  considerably  decreasing  the  image  quality, 
which  would  benefit  the  field  of  mammography  by  delivering  a  more  effective  balance  between  image 
quality  and  radiation  dose.  The  goal  of  this  study  was  therefore  to  provide  an  initial  investigation  of  the 
effects  of  a  range  of  beam  hardening  levels  on  the  DQE  and  the  average  glandular  dose. 
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Fig.  1.  Illustration  of  the  x-ray  configuration  utilized  in  this  study. 

3.  Materials  and  methods 

3.1.  Beam  hardening 

Beam  hardening  was  provided  through  the  use  of  the  commonly  used  mammography  filters  Rhodium 
(Rh)  and  Molybdenum  (Mo)  [3-5].  The  typical  filter  thicknesses  utilized  are  25  pm  of  Rh  and  30  pm  of 
Mo  [5].  This  study  thus  evaluated  four  filter  thicknesses  to  provide  a  thorough  comparison:  25  /tm  and 
50  pm  of  Rh,  and  30  //in  and  60  /im  of  Mo.  In  addition,  the  investigation  included  a  fifth  mode  without 
added  filtration  to  provide  a  baseline  for  the  comparison. 

3.2.  Inline  x-ray  imaging  system  and  measurement  components 

The  inline  x-ray  imaging  system  prototype  utilized  in  this  study  consists  of  the  imaging  and  measure¬ 
ment  components  mounted  along  an  optical  rail,  which  allows  the  use  of  a  precise  laser  alignment  process 
detailed  in  a  previous  study  [34].  This  also  provides  the  ability  to  adjust  the  x-ray  source-to-object  (Rl) 
and  object-to-detector  (R2)  distances  to  facilitate  the  x-ray  configuration  utilized  in  this  study,  which  is 
given  in  Fig.  1.  In  this  study,  Rl  and  R2  were  both  91.44  cm,  which  produced  a  source-to-detector  dis¬ 
tance  (R1+R2)  of  182.88  cm.  A  longer  distance  than  typical  clinical  mammography  values  was  utilized 
to  target  future  applications,  as  this  configuration  has  been  used  for  preliminary  studies  in  inline  phase 
contrast  x-ray  imaging  indicating  the  potential  to  reduce  the  dose  in  mammography  without  negatively 
affecting  the  detection  capability  [7,35-39].  In  addition,  the  air  gap  in  this  configuration  provides  a 
comparable  amount  of  scatter  rejection  and  resultant  image  quality  improvement  as  the  grid  used  in 
clinical  systems  [6,40,41].  The  object  was  simulated  by  a  BR12  phantom  (Model  014A,  CIRS,  Norfolk, 
Virginia,  USA)  with  a  thickness  of  5  cm.  BR12  phantoms  were  developed  to  simulate  a  50%  glandular, 
50%  adipose  tissue  composition  [42],  and  thus  are  used  extensively  for  dose  estimation  purposes  in 
clinical  and  research  applications  [2].  For  comparison  purposes,  a  second  set  of  measurements  were 
performed  without  the  object,  in  an  effort  to  evaluate  the  effect  of  the  object  on  the  measurements.  The 
prototype  system  employs  an  x-ray  tube  with  a  Mo  anode  and  a  Be  output  window  (UltraBright  Mi¬ 
crofocus  Source,  Oxford  Instruments,  Scotts  Valley,  California,  USA).  The  tube  operates  between  x-ray 
energies  of  20  and  60  kV  with  an  output  power  range  of  10  to  60  W.  The  focal  spot  diameter  varies  with 
power  output,  thus  a  constant  power  output  of  20  W  was  maintained  in  this  study  to  ensure  a  constant 
focal  spot  diameter  of  20  pm,  which  corresponds  to  a  tube  current  of  0.33  milliamperes  (mA).  The  image 
detection  system  employed  in  this  study  is  a  computed  radiography  system  with  mammography  plates 
(Regius  190,  Konica  Minolta  Medical  Imaging,  Wayne,  New  Jersey,  USA),  which  provides  a  pixel  pitch 
of  43.75  pm.  Data  linearization  was  performed  through  a  process  presented  previously  [26,31,43]  to 
allow  quantitative  DQE  analysis  of  the  image  data.  The  entrance  exposure  level  was  measured  with  a 
calibrated  ionization  chamber  (10X9-180  ionization  chamber,  Model  9095  measurement  system,  Radcal 
Corporation,  Monrovia,  California,  USA).  Results  from  a  previous  study  [24]  indicated  that  utilizing 
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the  mean  of  numerous  measurements  reduces  the  error;  therefore  five  measurements  at  each  mode  were 
acquired.  The  x-ray  spectrum  was  measured  through  the  use  of  an  x-ray  spectrometer  with  a  3  x  3  x 
1  mm3  CdTe  detector  (Amptek  Incorporated,  Bedford,  Massachusetts,  USA),  using  a  pair  of  tungsten 
collimators  provided  by  the  manufacturer. 

3.3.  Imaging  parameters 

An  x-ray  energy  of  60  kV  was  utilized  for  the  measurements  acquired  in  this  preliminary  study. 
A  higher  x-ray  energy  than  current  clinical  values  was  utilized  in  an  effort  to  target  future  clinical 
applications  in  phase  contrast  imaging.  Initial  investigations  in  phase  contrast  imaging  have  indicated 
the  potential  of  increasing  the  x-ray  energy  to  reduce  the  dose  without  impacting  the  image  quality  [6, 
39,44-46] .  The  ability  to  apply  the  technique  of  beam  hardening  to  high  energy  phase  contrast  imaging 
could  provide  an  additional  dose  reduction,  which  would  further  increase  the  benefits  to  the  field  of 
mammography.  However,  beam  hardening  alone  could  potentially  be  utilized  in  a  current  clinical 
environment  for  an  immediate  dose  reduction;  therefore,  subsequent  studies  will  further  evaluate  the 
feasibility  of  beam  hardening  through  comparisons  with  a  range  of  diagnostic  energies,  including  current 
clinical  mammography  values. 

The  exposure  time  for  each  filter/object  combination  was  selected  to  provide  a  constant  detector 
entrance  exposure  of  10  mR,  in  an  effort  to  provide  a  comparison  based  on  similar  radiation  exposure 
for  image  acquisition.  Once  the  exposure  time  was  determined  for  each  filter/object  combination,  the 
corresponding  object  entrance  exposure  was  measured  for  calculation  of  the  average  glandular  dose. 
Due  to  the  constant  detector  entrance  exposure  for  all  combinations,  the  differences  among  the  DQE  and 
dose  values  was  a  result  of  only  the  factors  controlled  by  this  study. 

3.4.  Average  glandular  dose  measurements 

As  demonstrated  in  Eq.  (1),  the  average  glandular  dose  Dg  is  based  on  the  object  entrance  exposure 
(. Xese )  and  the  average  glandular  dose  coefficient  (Dgjy).  The  object  entrance  exposure  measurements 
were  detailed  previously,  and  the  Dgiy  values  for  each  mode  were  estimated  with  Monte  Carlo  simulations 
through  a  process  detailed  in  previous  studies  [14,17,18].  The  simulations  assume  the  presence  of  an 
object  with  50%  adipose  and  50%  glandular  tissue  composition  in  the  path  of  the  x-ray  beam.  Therefore, 
the  Dgiy  values  were  calculated  only  for  the  spectra  acquired  when  the  object  is  utilized. 

3.5.  Detective  quantum  efficiency  measurements 

The  DQE  calculation  given  in  Eq.  (2)  consists  of  measurements  of  the  MTF,  NPS,  5(0)  and  q.  First, 
the  MTF  was  calculated  through  acquiring  one  image  of  a  10  pm  wide  slit  camera  (Nuclear  Associates, 
Carle  Place,  NY,  USA)  for  each  filter/object  combination.  The  MTF  of  each  image  was  determined 
according  to  previous  methods  [25,47-50]  by  normalizing  the  absolute  value  of  the  line  spread  function 
(LSF),  which  is  the  response  of  the  imaging  system  to  a  line  stimulus  [5].  Next,  the  noise  power  spectrum 
was  measured  through  acquisition  of  20  images  for  each  mode  in  the  absence  of  any  additional  objects 
in  the  path  of  the  x-ray  beam,  which  is  referred  to  as  a  noise-only  image.  Due  to  the  stochastic  nature  of 
noise  in  x-ray  images,  the  experimental  calculation  of  the  two-dimensional  (2D)  NPS  typically  involves 
separating  the  noise-only  image  into  numerous  smaller  regions  and  averaging  the  NPS  values  calculated 
from  each  region  [22,27,29,31,43,51],  The  one-dimensional  (ID)  NPS  is  then  determined  through  the 
use  of  a  slice  consisting  of  four  data  lines  parallel  to  and  immediately  adjacent  to  the  axes  [22,43], 
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Table  1 

Average  glandular  dose  calculations  for  an  x-ray  energy  of  60  kV 
and  20  W  with  a  5  cm  BR12  phantom,  assuming  50%  adipose,  50% 
glandular  tissue  composition.  The  exposure  time  was  selected  for 
each  filter  to  provide  a  constant  detector  entrance  exposure  of  10  mR. 
The  results  indicate  a  dose  reduction  for  all  filters 


Filtration 

Object 
exposure  (R) 

DgN 

(mrad/R) 

Dg 

(mrad) 

Dose 

reduction  (%) 

None 

0.81 

333.4 

268.69 

25  /im  Rh 

0.61 

343.8 

210.02 

22 

50  fim  Rh 

0.56 

366.6 

203.82 

24 

30  /im  Mo 

0.68 

337.0 

229.17 

15 

60  /im  Mo 

0.63 

356.0 

222.86 

17 

The  large  area  signal  5(0)  was  determined  by  calculating  the  mean  pixel  intensity  of  the  20  noise-only 
images  within  the  region  of  interest  utilized  for  calculation  of  the  noise  power  spectrum.  The  value  of  q 
is  determined  experimentally  through  multiplying  the  detector  entrance  exposure  by  the  photon  fluence 
per  unit  exposure,  which  can  be  calculated  from  a  measured  x-ray  spectrum  [24,43,52], 


4.  Results  and  discussion 

4. 1.  Average  glandular  dose 

Comparisons  of  the  average  glandular  dose  calculated  for  each  filter  with  the  5  cm  BR12  object  arc 
given  in  Table  1,  which  provides  the  measured  object  entrance  exposure  and  the  corresponding  range  of 
DgN  values  calculated  from  the  Monte  Carlo  estimations.  The  last  column  presents  the  percentage  of 
dose  reduction  provided  by  each  filter  as  compared  to  the  dose  value  obtained  without  filtration.  The 
results  demonstrate  a  notable  dose  reduction  of  15%  to  24%  for  the  range  of  beam  hardening  levels.  The 
filter  corresponding  to  the  maximum  dose  reduction  was  the  50  /mi  Rh  filter.  It  is  also  interesting  to 
note  that  the  filter  providing  the  second  largest  dose  reduction  was  the  25  /mi  Rh  filter,  which  provided 
a  lower  dose  than  both  Mo  filters. 

4.2.  Detective  quantum  efficiency 

Comparisons  of  the  x-ray  spectra  acquired  for  each  beam  hardening  level  arc  given  in  Fig.  2  (a)  and 
(b),  which  correspond  to  the  images  acquired  without  an  object  and  with  the  5  cm  BR12,  respectively. 
The  display  range  in  both  figures  has  been  adjusted  to  display  a  subset  of  the  full  x-ray  spectra  range, 
in  an  effort  to  more  effectively  demonstrate  the  differences  between  the  beam  hardening  levels  in  the 
region  of  interest  around  the  peaks.  The  behavior  in  the  low  and  high  energy  values  not  depicted  exhibits 
a  decreasing  Rend  to  0  in  both  directions.  The  comparison  without  the  object  demonstrates  the  effects 
of  beam  hardening  in  removing  low  energy  photons  from  the  x-ray  beam,  while  the  comparison  with  the 
object  illustrates  the  further  hardening  effects  of  the  object  on  the  x-ray  beam. 

Comparisons  of  the  MTF,  NPS,  and  NEQ  values  calculated  for  each  beam  hardening  level  are  given 
in  Figs  3  through  5,  respectively,  in  which  (a)  and  (b)  correspond  to  the  images  acquired  without  an 
object  and  with  the  5  cm  BR12,  respectively.  In  each  of  the  comparisons  without  the  object,  the  curve 
calculated  without  added  beam  hardening  is  noticeably  different  from  the  curves  corresponding  to  the 
filters.  However,  when  the  5  cm  BR12  is  utilized  to  simulate  the  human  breast,  the  separation  among 
the  curves  is  much  smaller,  which  is  an  encouraging  indication  that  the  use  of  filtration  in  a  clinical 
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Fig.  2.  Comparison  of  the  x-ray  spectra  corresponding  to  each  beam  hardening  level.  The  spectra  were  acquired  at  60  kV,  20  W 
(a)  without  an  object  in  the  path  of  the  x-ray  beam,  and  (b)  with  the  5  cm  BR12  object  in  the  path. 


Fig.  3.  Comparison  of  the  MTF  curves  corresponding  to  the  beam  hardening  levels.  The  curves  were  calculated  from  slit 
camera  images  acquired  at  60  kV,  20  W  (a)  without  an  object  in  the  path  of  the  x-ray  beam,  and  (b)  with  the  5  cm  BR12  object 
in  the  path. 

environment  may  not  negatively  affect  the  imaging  response  of  the  system.  For  example,  in  Fig.  3(a), 
the  shape  of  the  MTF  is  affected  by  each  level  of  added  filtration,  which  can  probably  be  attributed  to  the 
sensitivity  of  the  x-rays  to  the  filtration  without  an  object  in  the  path  of  the  x-ray  beam.  When  the  object 
is  utilized,  the  MTF  for  each  beam  hardening  level  is  degraded  in  comparison  to  the  corresponding  MTF 
without  the  object,  which  indicates  the  negative  effect  of  the  object.  This  is  intuitive,  due  to  the  amount 
of  attenuation  experienced  within  the  object.  In  this  comparison,  the  minimal  difference  between  the 
MTF  without  filtration  and  the  filtered  MTFs  indicates  that  beam  hardening  provides  a  limited  amount 
of  additional  degradation  to  the  imaging  response  of  the  system. 
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Fig.  4.  Comparison  of  the  NPS  curves  corresponding  to  the  beam  hardening  levels.  The  curves  were  calculated  from  noise-only 
images  acquired  at  60  kV,  20  W  (a)  without  an  object  in  the  path  of  the  x-ray  beam,  and  (b)  with  the  5  cm  BR12  object  in  the 
path. 


Fig.  5.  Comparison  of  the  NEQ  curves  corresponding  to  the  beam  hardening  levels.  The  curves  were  calculated  from  the  MTF, 
NPS  and  5(0),  (a)  without  an  object  in  the  path  of  the  x-ray  beam,  and  (b)  with  the  5  cm  BR12  object  in  the  path. 


Next,  comparisons  of  the  q  values  corresponding  to  the  beam  hardening  levels  are  given  in  Tables  2 
and  3,  which  correspond  to  the  measurements  acquired  without  an  object  and  with  the  5  cm  BR12, 
respectively,  in  the  path  of  the  x-ray  beam.  The  results  correlate  with  the  I)ny-  coefficients,  through 
exhibiting  an  increasing  trend  from  the  minimum  value  without  added  beam  hardening  to  the  maximum 
value  with  the  50  fim  Rh  filter,  which  also  produced  the  largest  D9n  value.  This  follows  intuition,  as 
the  largest  q  value  indicates  the  highest  concentration  of  photons  in  the  same  surface  area,  which  would 
also  correspond  to  the  largest  dose  coefficient. 

Finally,  comparisons  of  the  DQE  values  for  the  beam  hardening  levels  without  an  object  and  with 
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Table  2 

Comparison  of  the  photon  fluence 
(q)  values  for  the  beam  hardening 
levels,  corresponding  to  the  mea¬ 
surements  acquired  without  an  ob¬ 
ject  in  the  path  of  the  x-ray  beam 


Filtration 

Photon  fluence 
(q)  (1/mm2) 

None 

3.134  x  10s 

25  fim  Rh 

4.314  x  10s 

50  fim  Rh 

5.089  x  10s 

30  pi  Mo 

4.024  x  10s 

60  pi  Mo 

4.498  x  10s 

Table  3 

Comparison  of  the  photon  fluence 

(q)  values  for  the  beam  hardening 

levels,  corresponding  to  the  mea¬ 

surements  acquired  with  the  5  cm 

BR12  object  in  the  path  of  the  x- 

ray  beam 

Filtration 

Photon  fluence 

(q)  (1/mm2) 

None 

7.087  x  10s 

25  pi  Rh 

7.636  x  10s 

50  pi  Rh 

8.288  x  10s 

30  pi  Mo 

7.618  x  10s 

60  pi  Mo 

8.243  x  10s 

(a)  (b) 


Fig.  6.  Comparison  of  the  DQE  curves  corresponding  to  the  range  of  beam  hardening  levels,  (a)  without  an  object  in  the  path 
of  the  x-ray  beam,  and  (b)  with  the  5  cm  BR12  object  in  the  path. 

the  5  cm  BR12  are  provided  in  Fig.  6  (a)  and  (b),  respectively.  The  results  also  exhibit  much  closer 
correlation  between  the  curves  acquired  with  the  5  cm  BR12  than  the  curves  generated  without  the  object. 
Furthermore,  within  the  DQE  curves  corresponding  to  the  object,  once  again  only  a  small  divergence  is 
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distinguishable  between  the  curve  acquired  without  beam  hardening  and  the  curves  corresponding  to  the 
range  of  beam  hardening  levels. 


5.  Conclusion 

For  the  specific  experimental  conditions  utilized  in  this  preliminary  study,  the  results  of  the  dose  and 
DQE  comparisons  are  an  initial  indication  that  the  use  of  beam  hardening  holds  the  potential  to  reduce 
the  dose  without  decreasing  the  performance  of  the  system.  Subsequent  studies  will  further  evaluate 
the  feasibility  of  beam  hardening  through  a  thorough  investigation  of  any  differences  in  subject  contrast 
resulting  from  beam  hardening,  which  could  affect  the  image  quality  provided  by  the  system.  The 
ability  of  beam  hardening  to  more  effectively  balance  the  tradeoff  between  dose  and  image  quality  would 
benefit  the  fields  of  diagnostic  x-ray  imaging,  especially  mammography.  Future  studies  will  also  perform 
comparisons  encompassing  a  range  of  diagnostic  energies,  including  current  clinical  mammography 
values. 
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Purpose:  The  high  breast  density  is  one  of  the  biggest  risk  factors  for  breast  cancer.  Identifying 
patient  having  persistent  high  breast  density  is  important  for  breast  cancer  screening  and  prevention. 

In  this  work  the  authors  propose  for  the  first  time  an  x-ray  phase-shifts-based  method  of  breast  density 
measurement. 

Methods:  When  x  ray  traverses  the  breast,  x  ray  gets  not  only  its  intensity  attenuated  but  also  its 
phase  shifted.  Studying  the  x-ray  phase-shifts  generated  by  the  breast  tissues,  we  derived  a  general 
formula  for  determining  the  volumetric  breast  density  from  the  breast  phase  map.  The  volumetric 
breast  density  is  reconstructed  by  retrieving  the  breast  phase  map  from  just  a  single  phase-sensitive 
projection  of  the  breast,  through  the  use  of  an  innovative  phase  retrieval  method  based  on  the  phase- 
attenuation  duality.  In  order  to  numerically  validate  this  phase-shifts-based  method  for  measuring 
the  volumetric  breast  density,  the  authors  performed  computer  simulations  with  a  digitally  simulated 
anthropomorphic  breast  phantom. 

Results:  Using  the  proposed  phase-shifts-based  method,  we  reconstructed  the  breast  phantom’s  vol¬ 
umetric  breast  density,  which  differs  from  the  phantom’s  intrinsic  breast  density  by  only  0.06%. 

In  the  presence  of  noises  in  the  projection  image,  the  reconstructed  volumetric  breast  density  dif¬ 
fers  from  the  phantom’s  intrinsic  breast  density  by  only  1.79%  for  a  projection  signal-to-noise-ratio 
(SNR)  of  34.  The  error  in  reconstructed  breast  density  is  further  reduced  to  1.61%  and  1.55% 
for  SNR  =  68  and  SNR  =  134,  respectively,  achieving  good  accuracies  in  the  breast  density 
determination. 

Conclusions:  The  authors  proposed  an  x-ray  phase-shifts-based  method  of  measuring  the  volumetric 
breast  density.  The  simulation  results  numerically  validated  the  proposed  method  as  a  novel  method 
of  breast  density  measurement  with  good  accuracies.  ©  2012  American  Association  of  Physicists  in 
Medicine,  [http://dx.doi.org/10. 1 1 18/1 .4729838] 

Key  words:  X-ray  phase  contrast,  phase  retrieval,  breast  density 


I.  INTRODUCTION 

As  is  well  known,  breast  cancer  is  the  most  common  cancer 
among  American  women  and  is  the  leading  cause  of  cancer 
related  death  among  women  between  35  and  50  yr  of  age 
in  the  United  States.  In  recent  years  scientific  studies  have 
found  that  the  breast  cancer  risk  is  directly  linked  to  the  breast 
density,  the  relative  amount  of  the  dense  fibroglandular  tissue 
versus  fatty  adipose  tissue,  namely,  the  fraction  of  the  vol¬ 
ume  of  fibroglandular  tissue  in  the  breast.  As  the  fibroglan¬ 
dular  tissues  appear  “hyperdense”  in  mammograms  because 
of  its  relatively  high  attenuation  coefficients,  studies  found 
that  high  breast  density  inferred  from  mammograms  is  as¬ 
sociated  with  twofold  to  fourfold  increase  in  risk  for  breast 
cancer. 1-3  Although  the  magnitude  of  that  risk  is  still  under 
some  debate,  the  scientific  literature  holds  that  high  breast 
density  is  one  of  the  biggest  risk  factors  for  breast  cancer.2 
Identifying  patient  having  persistent  breast  density  may  ben¬ 
efit  these  patients  from  a  more  frequent  breast  cancer  screen¬ 


ing  for  early  detection  of  breast  cancer  or  from  some  other 
preventive  measure  such  as  the  chemoprevention.  The  pub¬ 
lic  awareness  of  breast  density  as  a  high  risk  factor  for  breast 
cancer  has  recently  been  pushing  new  legislature  and  regu¬ 
latory  efforts  for  providing  patients  the  breast  density  infor¬ 
mation  from  their  breast  imaging  exams.  Therefore,  in  recent 
years  there  is  strong  interest  in  developing  methods  to  mea¬ 
sure  volumetric  breast  density  (VBD).  The  most  straightfor¬ 
ward  method  of  measuring  volumetric  breast  density  would 
be  from  the  three-dimensional  x-ray  breast  image  techniques 
such  as  breast  computed  tomography  (CT),  which  is  still  in 
its  early  period  of  development.  Currently,  it  is  most  prac¬ 
tical  to  determine  volumetric  breast  density  from  the  two- 
dimensional  imaging  techniques  such  as  mammography. 4-1(1 
The  mammography-associated  methods  are  all  based  on  the 
small  differences  in  tissue  attenuation  of  the  dense  breast  fi¬ 
broglandular  tissue  and  fatty  adipose  tissue.  In  some  methods, 
breast  density  is  measured  by  the  categorical  scores  based  on 
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the  tissue  radiographic  opacity  on  a  mammogram  acquired 
with  a  screen/film  based  or  digital  system.4,9' 10  However, 
such  methods  are  limited  by  lack  of  consensual  reference 
standards  and  are  often  flawed  by  the  subjectivity.  The  im¬ 
proved  methods  reported  in  the  literature  all  need  to  con¬ 
duct  elaborated  phantom-based  system  calibrations  for  breast 
density  measurement.  Some  of  these  methods  require  simulta¬ 
neous  exposure  of  a  calibration  phantom  alongside  the  breast 
for  the  system  calibrations.7, 8,10  These  laborious  calibrations 
are  indeed  required,  because  mammography  exams  are  con¬ 
ducted  with  different  tube  targets/filters  and  voltages,  and  the 
linear  attenuation  coefficients  of  breast  tissues  vary  in  compli¬ 
cated  ways  as  the  x-ray  energy  changes.  Therefore,  in  search 
for  a  more  effective  method  of  breast  density  measurement, 
it  would  be  desirable  to  develop  a  method  to  measure  breast 
density  directly  rather  than  measuring  tissue  attenuation  as  a 
surrogate  of  the  breast  density. 

In  this  work  we  propose  to  determine  the  breast  density 
by  direct  probing  the  projected  electron  densities  of  breast 
fibroglandular  tissue  and  adipose  tissue  through  the  use  of 
the  retrieved  phase  map  of  a  breast.  Note  that  the  elec¬ 
tron  density  of  a  tissue  is  the  intrinsic  property  of  the  tis¬ 
sue  and  it  is  independent  of  x-ray  photon  energies  employed. 
In  Secs.  II-IV,  we  lay  down  a  framework  of  a  phase-shifts- 
based  method  of  determining  the  volumetric  breast  density. 
We  show  how  to  retrieve  the  breast  phase  map  from  a  single 
phase-sensitive  projection  of  the  breast.  In  order  to  validate 
this  novel  breast  density  measurement  method,  we  discuss  the 
volumetric  breast  density  values  reconstructed  from  computer 
simulations  with  a  digitally  simulated  anthropomorphic  breast 
phantom. 

II.  METHODS 

Our  method  of  breast  density  measurement  is  based  on  an 
x-ray  phase-sensitive  imaging  technique.  As  x  rays  traverse 
the  breast,  x  rays  undergo  phase-shifts  in  addition  to  the  inten¬ 
sity  attenuation.  The  amount  of  x-ray  phase-shift  ip  generated 
by  tissues  is  equal  to  <p  =  — ( X)refmype(s)ds  =  —  {X)repe^p, 
where  re  =  2.818  x  1  0  15  m  is  the  classical  electron  radius 
and  X  denotes  the  x-ray  wavelength.1 1-1  s  In  this  paper  the 
angle  bracket  ( •  >  denotes  the  spectral  average  over  the  de¬ 
tected  polychromatic  x-ray  spectrum,  so  (X)  denotes  the  av¬ 
erage  x-ray  wavelength.  Here,  pe  denotes  the  tissue  electron 
density  and  pep  =  fraype(s)ds  is  the  projected  electron  den¬ 
sity  summed  over  the  x-ray  path.  Therefore,  once  we  have  re¬ 
trieved  the  phase  map  q>  (?)  of  a  breast  through  its  projection 
image,  where  ?  denotes  the  position  on  the  phase  map,  we  are 
able  to  find  the  projected  electron  density  map  pep  (?)  of  the 
breast  as  pe,P(r )  =  (— 1/(A .)re)<p(r).  Using  the  pep(?)-inap  of 
a  breast,  we  can  determine  the  volumetric  breast  density  of 
the  breast  as  follows.  For  the  purpose  of  breast  density  deter¬ 
mination,  we  may  assume  that  there  are  only  two  types  of  tis¬ 
sues  in  breast:  the  fibroglandular  tissue  with  an  electron  den¬ 
sity  pe  fg  and  the  adipose  tissue  with  an  electron  density  pe  a d. 
Note  that  the  values  of  pe^g  and  pe  ad  are  available  from  the 
measured  data  in  the  literature,  pe$g  =  3.448  x  1023/cm3  and 
Pe.ad  =  3.108  x  1023/cm3.16  Let  m  (?)  be  the  fibroglandular- 


tissue  fraction  along  a  projection-ray  detected  at  position  ?, 
then  we  have 


Pe,p(r)  =  [m(?)  •  peSg  +  (1  -  m(r))  ■  peM\  ■  T(r), 


(1) 


where  T  (?)  denotes  the  ray  path-length  in  the  breast.  In  mam¬ 
mography,  the  breast  is  always  kept  forcefully  and  uniformly 
compressed  for  spreading  breast  tissues  and  reducing  the  scat¬ 
ters.  The  compressed  breast  thickness  is  always  automati¬ 
cally  measured  and  indicated  by  the  imaging  system.  As  long 
as  the  source-breast  distance  is  significantly  larger  than  the 
breast  lateral  size,  the  values  of  T  (?)  are  approximately  equal 
to  the  compressed  thickness  Tc  of  the  breast,  except  for  the 
breast  anterior  edge,  which  is  not  in  contact  with  the  com¬ 
pression  paddle  and  constitutes  only  a  very  small  portion  of 
the  breast  volume.  With  this  approximation  we  found  that  the 
fibroglandular-tissue  fraction  along  the  ray-?  as  given  by 

~(p(r)/{X)re  -  pe,a&Tc 


m(r) 


1  1  Pe,  ad)Tc 


(2) 


While  m  (?)  specifies  the  fibroglandular-tissue  fraction  along 
the  ray-?,  the  breast  density,  as  a  risk  factor  defined  in  the 
literature,  should  be  quantified  the  fibroglandular-tissue  frac¬ 
tions  averaged  over  the  whole  breast.  Hence,  a  map  of  these 
m  (?)-values,  which  we  call  the  m  (?)-image  of  a  breast,  is 
the  map  of  the  breast’s  fibroglandular-tissue  fractions.  The 
VBD  can  then  be  computed  from  the  m  (?)-image  as  VBD 
=  Abreast  m(r)/Np,  where  the  sum  is  taken  over  the  projected 
breast  image  and  Np  is  the  total  pixel  number  in  the  projected 
breast.  Hence,  the  volumetric  breast  density  can  be  found 
from  the  breast  phase  map  as 

Ebreast  (~<P(r)/ {X)re  -  pe,adTc) 


VBD  = 


NpTc(pe,f g  Pe,  ad) 


(3) 


In  order  to  implement  a  phase-shifts-based  method  of  mea¬ 
suring  the  breast  density,  it  calls  for  a  low  radiation  dose  yet 
phase-sensitive  technique  to  recover  the  phase-shifts  gener¬ 
ated  by  the  breast,  as  the  average  glandular  doses  involved 
should  be  closely  controlled  and  monitored  as  is  in  current 
mammography.  Although  the  differences  in  x-ray  phase-shifts 
between  different  tissues  are  about  1000  times  greater  than 
their  differences  in  the  attenuation,  current  technology  is  un¬ 
able  to  directly  measure  the  phase-shifts  because  of  the  ex¬ 
tremely  high  frequencies  of  x  rays.1 1-13  Breast  phase-shifts 
have  to  be  retrieved  from  the  phase-sensitive  images  of  the 
breast.  The  proposed  Eq.  (3)-based  method  of  breast  density 
measurement  can  be  applied  to  any  of  phase-sensitive  imag¬ 
ing  techniques  such  as  the  inline  or  the  grating-based  tech¬ 
niques.  Among  these  phase-sensitive  techniques,  the  inline 
phase-sensitive  imaging  is  the  simplest  to  be  implemented, 
since  the  inline  technique  does  not  need  any  x-ray  optics  such 
as  the  crystal  analyzers  or  high  line-density  transmission  grat¬ 
ings.  In  this  work,  we  study  the  inline  phase-sensitive  imaging 
only.  The  setting  for  the  inline  technique  is  similar  to  that  of 
conventional  mammography,  provided  one  uses  a  small  focal 
spot  and  a  sufficiently  large  object-detector  distance.12,14'15 
In  the  inline  phase-sensitive  imaging,  as  x  rays  traverse  the 
breast,  x  rays  undergo  phase-shifts  and  intensity-attenuation, 
and  then  diffract  freely  over  a  distance  upon  detection.  The 
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resulted  image  exhibits  not  only  the  tissue  attenuation  con¬ 
trast,  but  also  the  tissue  phase  contrast,  which  manifests  as 
the  dark-bright  fringes  at  tissues’  boundaries  and  interfaces  in 
the  image.  11-15  In  order  to  measure  the  breast  density  by  us¬ 
ing  the  phase  map  of  a  breast,  one  needs  to  retrieve  the  breast 
phase  map  from  its  phase-sensitive  projection  images. 

The  phase  retrieval  is  based  on  x-ray  propagation  equa¬ 
tion,  which  reveals  how  the  phase-shifts  are  encoded  in  the 
image  intensity  variations.  A  convenient  form  of  the  propaga¬ 
tion  equation  is  the  so-called  transport  of  intensity  equation 
(TIE):17-18 


^  =  W- 


A2a(rD/M) 


R2W 

'  2nM 


V  •  {A20(rD/M)V<p(rD/Mj) 


(4) 


where  I  ( ro )  is  the  detected  x-ray  intensity,  /,,,  is  the  entrance 
x-ray  intensity,  A2(r)  is  the  attenuation  image  of  the  imaged 
object,  say  a  breast,  and  <p(r)  is  the  x-ray  phase-shift  map  of 
the  breast.  In  Eq.  (4)  /A  is  the  object-detector  distance,  M  is 
the  magnification  factor  employed  in  the  projection,  and  r o 
denotes  the  position  in  the  detector  plane.  The  operator  V  de¬ 
notes  the  two-dimensional  transverse  gradient  differential  op¬ 
erator.  From  Eq.  (4)  it  is  clear  that  the  detected  x-ray  intensity 
I  (ro)  is  determined  not  only  by  the  conventional  attenuation- 
based  contrast  A 2  (r)  but  also  by  the  transverse  Laplacian  and 
gradient  differentials  of  the  phase-shift  (p  (r).  These  differen¬ 
tials  represent  the  edge-enhancements  exhibit  in  the  phase- 
sensitive  projection  images.  Since  the  phase  contrast  and  at¬ 
tenuation  contrast  are  mixed  together  in  a  phase-sensitive  pro¬ 
jection,  there  are  two  unknown  variables  for  a  given  projec¬ 
tion  intensity  in  Eq.  (4).  Hence,  the  TIE-based  phase  retrieval 
method  requires  multiple  projections  (at  least  two  projections) 
acquired  with  varying  object-detector  distances  for  retrieving 
the  phase-shift  map  (p(r )  of  a  breast. 17-19  This  requirement 
of  multiple  image-acquisitions  for  the  TIE-based  phase  re¬ 
trievals  is  obviously  cumbersome  in  implementation,  and  the 
multiple  exposures  make  the  radiation  doses  to  breast  multi¬ 
plied.  As  if  these  were  not  enough  challenges,  the  TIE-based 
phase  retrieval  method  is  indeed  unstable  against  the  noise 
in  projection  images.  While  one  might  want  to  increase  the 
radiation  doses  employed  in  the  projections  for  retrieving  a 
quality  phase  map  of  a  breast,  the  glandular  dose  to  a  breast 
should  be  kept  as  low  as  reasonably  achievable  because  the 
breast  is  a  radiation  sensitive  organ.  Considering  these  disad¬ 
vantages  with  the  TIE-based  phase  retrieval  method,  we  have 
to  employ  a  more  effective  phase  retrieval  method  of  measur¬ 
ing  breast  density  through  the  use  of  Eqs.  (2)  and  (3). 

We  noted  that  the  breast  fibroglandular  tissues  and  adipose 
tissues  are  consisting  of  dominantly  elements  with  Z  <  10. 
If  we  employ  high  energy  x  rays  of  60  keV  or  higher  for 
imaging,  the  x-ray-breast  tissue  interactions  are  dominated  by 
the  x-ray  Compton  scattering  from  freelike  atomic  electrons, 
while  the  x-ray  photoelectric  absorption  and  coherent  scatter¬ 
ing  are  all  diminished  and  negligible.21’  In  this  situation,  both 
the  tissue  attenuation  and  phase-shift  are  all  determined  by 


tissues’  electron  density  distributions  such  that 
A2(r)  =  exp[  ctkn  (E)  pe,P(r )], 

<p(r)  =  -^repe<p(r),  (5) 


where  cr  knCE)  denotes  the  Klein-Nishina  total  cross  section 
for  Compton  scattering  and  E  is  the  photon  energy  The  Klein- 
Nishina  total  cross  section  changes  slowly  with  photon  energy 


ctkn  (E)  =  2nr2 


TjL 


2(1  +  77)  1 

iogd  +  21,) 

1+2  rj  r] 


1  1+3)7 

-l0g(l  +  2„)-^-^ 


(6) 


where  i]  =  E/511  keV.  In  a  recent  work  we  explored  the 
relationship  between  the  phase- shift  and  attenuation  defined 
by  Eq.  (5)  and  called  this  as  the  phase-attenuation  duality 
(PAD).21  We  found  that  when  the  phase-attenuation  duality 
holds,  the  x-ray  propagation  equation  gets  simplified  and  the 
phase  map  can  be  retrieved  from  just  a  single  phase-sensitive 
projection  I  (ro)  of  the  imaged  object:21'22 

<P(?)  =  ln{[l  -  (tf/a^Rir^/hrM)]-1 
(°kn) 

■(M2I(7o)/hn)}  (7) 


Here,  the  operator  V2  denotes  the  two-dimensional  transverse 
Laplacian  differential  operator  V2  =  (d2/dx2  +  d2/d y2)  and 
the  operator  D(V2)  =  [1  —  ({X2 /<Jkn)  ^2AV2 /27r M)]-1  is  a 
pseudodifferential  operator.  The  action  on  a  function  g(r)  of 
a  pseudodifferential  operator  such  as  V  (V2)  is  defined  as 

V(V2)g(T)  =  JJ  exp(/ 2n (s  -?)•/)•  D(-4 n2]2) 

■g(s)d2s-d2f,  (8) 

where  s  and  /  denote  the  integral  variables  in  the  transverse 
coordinate  space  and  frequency  space,  respectively.  We  call 
Eq.  (7)  the  PAD-based  phase  retrieval  formula.21  Recently, 
we  have  theoretically  and  experimentally  demonstrated  that 
the  PAD-based  phase  retrievals  are  robust  against  the  noise, 
in  a  stark  contrast  to  the  TIE-based  method  that  is  unstable 
in  the  presence  of  noise  in  images.22  Of  course,  in  order  for 
breast  tissues  to  fulfill  the  phase-attenuation  duality  condi¬ 
tion,  the  tube  voltages  should  be  set  to  as  high  as  120-150 
kVp,  grossly  different  from  the  tube  voltages  employed  in  cur¬ 
rent  mammography,  where  tube  voltages  used  are  lower  than 
40  kVp.  Although  the  attenuation-based  tissue  contrast  gets 
much  reduced  at  120-150  kVp,  the  phase-contrast  extracted 
from  the  phase  retrieval  will  compensate  for  the  attenuation- 
contrast  loss  associated  with  the  use  of  high  tube  voltages.22 
Putting  above  discussion  together,  we  propose  that  it  is  fea¬ 
sible  to  measure  the  projected  electron  densities  of  a  breast 
by  using  the  high-energy  x-ray  phase-sensitive  imaging  tech¬ 
nique.  One  can  retrieve  the  phase  map  of  a  breast  by  using 
only  a  single  phase-sensitive  projection  of  the  breast  and  per¬ 
forming  the  PAD-based  phase  retrieval  as  is  shown  in  Eq.  (7). 
From  the  retrieved  phase  map  of  the  breast  one  can  determine 
its  volumetric  breast  density  by  using  Eq.  (3). 
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In  order  to  numerically  validate  the  proposed  new  method 
of  breast  density  measurement,  we  constructed  a  digitally 
simulated  anthropomorphic  breast  phantom.  The  phantom  in¬ 
cludes  the  simulated  breast  ductal  system,  the  lobule  sys¬ 
tem,  and  three  masses  simulating  the  invasive  ductal  car¬ 
cinoma,  and  four  groups  of  small  calcification  spheres  (in 
sizes  of  0.54,  0.4,  0.32,  and  0.2  mm)  simulating  the  breast 
microcalcifications.  In  the  breast  phantom  each  voxel  was 
randomly  assigned  with  its  electron  density  pe  according  to 
a  random  Gaussian  process  with  a  zero  mean  and  unit  stan¬ 
dard  deviation.  To  simulate  the  breast  tissue  texture,  the  spa¬ 
tial  power  spectrum  of  the  random  electron  density  distribu¬ 
tion  was  further  modified  such  that  it  varies  as  |/|L35  for  a 
frequency  f .  The  resulted  phantom’s  intrinsic  by  construc¬ 
tion  volumetric  breast  density,  which  we  denote  by  VBD  |ruc, 
can  be  computed  from  the  assigned  voxel  electron  densities 
{pe( v),  v  =  1, . . . ,  N}  as 

N 

TV  Pe  Pe, ad 

VBDxrue  =  — ^ - ,  (9) 

Pe,  fg  Pe.ad 

where  N  is  the  total  number  of  the  phantom  voxels. 

In  order  to  numerically  validate  the  phase-based  method 
of  breast  density  measurement,  we  first  computed  the  breast 
projections  with  different  photon  energies.  In  the  simula¬ 
tions  we  assumed  a  detector  of  50  pm  pixels,  and  set  Pi 
=  R2  =  1  m.  The  source’s  focal  spot  was  assumed  as  a 
point  source.  This  is  an  adequate  approximation  as  long  as 
the  blur  caused  by  the  focal  spot  is  less  than  the  detector 
pixel  size.1’’  In  our  simulation  the  magnification  factor  M  = 
2  and  pixel  size  is  a  50  pm,  so  a  focal  spot  of  50  pm  or 
smaller  in  size  can  be  approximately  treated  as  a  pointlike 
source  in  the  simulation.  In  this  simulation  we  assumed  that 
the  source  was  with  a  tungsten  target  and  operated  at  140 
kVp.  We  assumed  that  the  detection  energy  spectrum  ranges 
from  54  keV  to  140  keV,  including  the  tungsten  characteris¬ 
tic  radiations.  Selecting  this  photon-energy  range  in  the  sim¬ 
ulation  was  to  satisfy  the  phase-attenuation  duality  condi¬ 
tion  with  good  approximation.  This  implied  that  the  source 
should  operate  at  high  voltages  and  with  heavy  filtrations 
in  the  implementation  of  the  proposed  method.  For  a  given 
photon  energy,  a  ray-tracing  algorithm  was  used  to  com¬ 
pute  the  ray  integrals  such  as  Jray  pe(r ,  s)ds  in  the  projection. 
The  x-ray  phase- shifts  and  attenuations  through  the  phan¬ 
tom  were  determined  as  <pph(^)  =  —  kre  /ray  pe(r ,  s)ds  and 
A2h(r)  =  exp[—  (TknCF1)  /ray  pe(r,  s)ds],  as  we  assumed  that 
the  phase-attenuation  duality  condition  is  satisfied  for  fibrog- 
landular  and  adipose  tissues  in  this  simulation  to  numerically 
validate  the  proposed  method.  Using  the  computed  cpph(r)  and 
Aph(r),  we  simulated  the  x-ray  Fresnel  diffraction  of  the  exit 
x  rays  propagating  from  the  breast  phantom  to  the  detector,23 
and  thereby  we  computed  the  detected  Fresnel  diffraction  in¬ 
tensity  /monof/D,  E)  for  each  individual  energy  E.  To  com¬ 
pute  the  phase-sensitive  projection  intensity  /poiy(?z>)  with  the 
polychromatic  x  rays,  we  summed  the  /lnono(f0,  E)  over  pho¬ 
ton  energies  in  the  detection  spectrum.  On  the  other  hand,  it  is 
important  to  examine  the  effects  of  the  noise  on  the  accuracies 


in  the  breast  density  measurement.  In  fact,  the  robustness  of 
the  phase  retrieval  method  is  critical  to  the  success  of  the  pro¬ 
posed  phase-based  method  of  breast  density  measurement,  as 
some  methods  in  the  literature  may  be  unstable  at  some  sub¬ 
stantial  noise  levels.22  To  simulate  the  noise  effects,  we  added 
the  Poisson  noise  with  different  noise-levels  into  the  simu¬ 
lated  phase-sensitive  projections  of  the  phantom,  following 
the  standard  noise  simulation  method  used  in  the  literature.24 
Employing  the  PAD-based  phase-retrieval  formula  (7),  we  re¬ 
trieved  the  phase  map  (p(r)  of  the  breast  phantom.  Using  the 
retrieved  (p(r)  and  Eq.  (3)  we  reconstructed  the  phantom’s  vol¬ 
umetric  breast  density  VBD. 

III.  RESULTS 

According  to  the  aforementioned  method,  we  constructed 
a  digitally  simulated  anthropomorphic  breast  phantom,  which 
simulates  a  compressed  breast  of  4.5  cm  thick,  20  cm  in  lat¬ 
eral,  and  10  cm  in  chest  wall-nipple  dimension.  Applying 
Eq.  (9)  to  the  randomly  assigned  electron  densities  of  this 
phantom,  we  found  the  phantom’s  intrinsic  by  construction 
volumetric  breast  density  VBDxrue  =  0.5025.  We  take  this 
VBDxrue  value  as  the  reference  to  evaluate  the  accuracy  of 
the  reconstructed  phantom  volumetric  breast  density. 

From  the  simulated  phase-sensitive  projection  image  of 
the  phantom  and  using  Eqs.  (2)  and  (7)  with  Tc  =  4.5  cm, 
we  computed  the  m(r)-image  of  the  breast  phantom,  the  map 
of  the  phantom’s  fibroglandular-tissue  fractions,  as  is  shown 
in  Fig.  1.  In  this  image,  the  shade  of  gray  of  a  pixel  repre¬ 
sents  the  computed  fraction  of  fibroglandular  tissues  along 
a  ray  passing  the  pixel.  A  lighter  shade  of  gray  represents 
a  larger  fraction  of  fibroglandular  tissues.  In  this  way,  the 
/«(r)-image  exhibited  good  breast  tissue  contrast  as  well 
(Fig.  1.)  Moreover,  this  m  (r)-image  provided  the  basis  for 
reconstructing  the  phantom’s  volumetric  breast  density.  Us¬ 
ing  the  VBD-formula  (3),  we  found  the  phantom’s  volumet¬ 
ric  breast  density  VBD  =  0.5028,  which  differs  from  its 
intrinsic  volumetric  breast  density  by  only  0.06%,  achiev¬ 
ing  an  excellent  accuracy.  Furthermore,  testing  the  breast 
density  determination  from  the  projection  images  with  dif¬ 
ferent  signal-to-noise-ratio  (SNR)  levels  associated  with 
different  mean  glandular  doses.  The  SNR,  which  is  simi¬ 
larly  defined  as  in  conventional  mammography,  is  defined 


Fig.  1 .  The  computed  m  (r)-image  of  the  breast  phantom,  namely,  the  map 
of  the  phantom’s  fibroglandular-tissue  fractions. 
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as  the  mean  intensity  over  a  ROI  drawn  in  the  center  of 
the  phantom  projection  image  over  the  standard  deviation  of 
the  intensity  values.  We  tested  the  cases  with  the  projection 
SNR  =  34,  64,  and  136,  which  are  corresponding  to  esti¬ 
mated  mean  glandular  doses  of  0.14,  0.57,  and  2.21  mGy, 
respectively.  We  found  that  the  reconstructed  VBD  values 
keep  achieving  good  accuracies,  in  spite  of  the  presence  of 
noise  in  the  simulated  projections.  With  the  projection  SNR 
=  34  we  found  the  reconstructed  VBD  =  0.5115,  which 
differs  from  the  intrinsic  volumetric  breast  density  VBDxme 
by  only  1.79%.  The  error  in  reconstructed  VBD  decreases 
with  increasing  projection  SNR  value.  For  example,  we  found 
the  reconstructed  VBD  =  0.5106  for  SNR  =  68,  which  differs 
from  VBDxme  by  only  1.61%.  The  error  in  VBD  is  further  re¬ 
duced  to  1.55%  for  SNR  =  134.  All  these  results  numerically 
validated  the  proposed  phase-shifts-based  method  as  a  novel 
breast  density  measurement  method  with  good  accuracies. 

IV.  DISCUSSION  AND  CONCLUSIONS 

Current  x-ray-based  methods  of  breast  density  measure¬ 
ment  are  all  based  on  small  attenuation  differences  between 
breast  tissues.  However,  as  the  linear  attenuation  coefficients 
of  breast  tissues  vary  in  a  complicated  way  as  the  employed 
x-ray  spectrum  changes,  all  these  methods  have  to  require 
elaborated  phantom-based  system  calibrations,  some  of  them 
include  simultaneous  exposures  of  the  phantom  alongside 
the  breast  for  system  calibrations.7  8, 10  The  proposed  phase- 
shifts-based  method  measures  the  tissue  projected  electron 
densities  rather  than  their  projected  attenuation  coefficients. 
Our  strategy  of  direct  probing  tissue  electron  densities  simpli¬ 
fies  the  system  calibration  for  breast  density  measurement.  In 
fact,  as  Eqs.  (3)  and  (7)  involve  the  spectral  averages  of  x-ray 
wavelengths  and  the  Klein-Nishina  cross  sections,  the  system 
parameters  that  our  method  requires  are  the  source  x-ray  spec¬ 
tra  and  detector  responses.  Another  key  idea  in  our  method  is 
to  use  the  PAD-based  phase  retrieval  technique  for  obtaining 
the  breast  phase  map.  As  breast  imaging  is  stringently  regu¬ 
lated  for  radiation  doses  involved,  there  are  always  substantial 
noises  in  breast  projection  images.  Some  common  phase  re¬ 
trieval  methods  in  the  literature  are  actually  unstable  in  pres¬ 
ence  of  noise,22  and  the  large  errors  in  the  retrieved  phase  map 
caused  by  the  retrieval  instability  would  spoil  the  breast  den¬ 
sity  measurement.  The  high  robustness  against  noise  of  the 
PAD-based  phase  retrieval  technique  is  crucial  to  our  method 
for  achieving  good  accuracies  in  breast  density  determination. 
Of  course,  many  future  works  are  required  for  implementing 
this  phase-based  method  of  breast  density  measurement.  In 
our  compressed  breast  model,  13.5%  of  breast  volume  is  in 
the  curved  periphery  zone  of  the  breast,  and  breast  thickness 
in  the  periphery  is  not  uniform  and  falls  off  quickly.  This  por¬ 
tion  of  tissues  was  not  included  in  the  VBD  computation  for 
avoiding  the  errors  in  thickness  estimation.  Therefore,  if  the 
periphery  has  a  quite  different  breast  density  than  the  rest 
of  breast,  this  approach  may  cause  errors  in  the  computed 
VBD.  Hence,  how  to  accurately  estimate  the  breast  thickness 
is  an  area  for  improvement,  and  it  is  also  an  active  research 
area  in  implementing  other  x-ray  imaging-based  methods  of 


breast  density  measurement.23  In  addition,  it  will  be  useful 
to  develop  a  correction  method  for  the  VBD  errors  in  cases 
where  an  imaging  setup  does  not  accurately  satisfy  the  phase- 
attenuation  condition. 

In  summary,  in  this  work  we  propose  for  the  first  time  an 
x-ray  phase-shifts-based  method  of  measuring  the  volumetric 
breast  density.  The  proposed  Eq.  (3)-based  method  of  breast 
density  measurement  can  be  applied  to  any  of  phase-sensitive 
imaging  techniques  such  as  the  inline  or  the  grating-based 
techniques.  This  phase-shifts-based  method  measures  the  tis¬ 
sue  projected  electron  densities  rather  than  their  projected 
attenuation  coefficients.  The  direct  probing  of  tissue  elec¬ 
tron  densities  in  our  method  simplifies  the  system  calibration 
requirement.  The  robust  PAD-based  phase  retrieval  method 
employed  in  the  inline  technique  ensures  the  accuracy  of 
the  reconstructed  volumetric  breast  densities.  Constructing  an 
anthropomorphic  digitally  simulated  breast  phantom,  we  per¬ 
formed  simulations  for  validating  this  novel  method  of  breast 
density  measurement.  Testing  with  different  projection  SNR- 
levels  in  the  simulations,  we  found  that  the  reconstructed  vol¬ 
umetric  breast  densities  achieved  good  accuracies  compared 
to  the  phantom’s  intrinsic  volumetric  breast  density.  All  these 
results  numerically  validated  the  proposed  phase-shifts-based 
method  as  a  novel  method  of  breast  density  measurement  with 
good  accuracies. 

ACKNOWLEDGMENTS 

This  research  was  supported  in  part  by  the  Department 
of  Defense  Breast  Cancer  Research  Program  under  Award 
No.  W81XWH-08-1-0613  and  the  NIH  under  Grant  No. 
R01CA142587.  H.  Liu  would  like  to  acknowledge  the  support 
of  Charles  and  Jean  Smith  Chair  endowment  fund  as  well. 

a) Author  to  whom  correspondence  should  be  addressed.  Electronic  mail: 
xwu@uabmc.edu 

*W.  J.  Nk.  "Risk  for  breast  cancer  development  determined  by  mammo- 
graphic  parenchymal  pattern,”  Cancer  37,  2486-2492  (1976). 

2N.  F.  Boyd,  L.  J.  Martin,  M.  J.  Yaffe,  and  S.  Minkin,  "Mammographic  den¬ 
sity  and  breast  cancer  risk:  Current  understanding  and  future  prospects,” 
Breast  Cancer  Res.  Treat.  13,  1-12  (201 1). 

3V.  A.  McCormack  and  I.  dos  Santos  Silva,  “Breast  density  and  parenchy¬ 
mal  patterns  as  markers  of  breast  cancer  risk:  A  meta-analysis,”  Cancer 
Epidemiol.  Biomarkers  Prev.  15,  1159-1169  (2006). 

4R.  P.  Highnam  and  J.  M.  Brady,  Mammographic  Image  Processing  (Kluwer 
Academic,  Boston,  1999). 

5R.  Highnam,  M.  Brady,  and  B.  Shepstone,  “A  representation  for  mammo¬ 
graphic  image  processing,”  Med.  Image  Anal.  1,  1-8  (1996). 

6J.  A.  Shephard.  K.  M.  Kerlikowske,  R.  Smith-Bindman,  H.  K.  Genant,  and 
S.  R.  Cummings.  "Measurement  of  breast  density  with  dual  X-ray  absorp¬ 
tiometry:  Feasibility,”  Radiology  223,  554-557  (2002). 

70.  Pawluczyk,  B.  J.  Augustine,  M.  J.  Yaffe,  D.  Rico,  J.  Yang,  G.  E.  Mawd- 
sley,  and  N.  F.  Boyd,  "A  volumetric  method  for  estimation  of  breast  density 
on  digitized  screen-film  mammograms,”  Med.  Phys.  30,  352—364  (2003). 
8S.  Van  Engeland,  P.  R.  Snoeren,  H.  Huisman,  C.  Boetes,  and  N.  Karsse- 
meijer,  "Volumetric  breast  density  estimation  from  full-field  digital  mam¬ 
mograms,”  IEEE  Trans.  Med.  Imaging  25,  273—282  (2006). 

9K.  E.  Martin,  M.  A.  Helvie,  C.  Zhou,  M.  A.  Roubidoux,  J.  E.  Bai¬ 
ley,  C.  Paramagul,  C.  E.  Blane,  K.  A.  Klein,  S.  S.  Sonnad,  and 
H.-P.  Chan,  "Mammographic  density  measured  with  quantitative  computer- 
aided  method:  Comparison  with  radiologists’  estimates  and  BIRADS  cate¬ 
gories,”  Radiology  240,  656-665  (2006). 

10M.  J.  Yaffe,  “Measurement  of  mammographic  density,”  Breast  Cancer  Res. 
Treat.  10,  1-10  (2008). 


Medical  Physics,  Vol.  39,  No.  7,  July  2012 


Wu,  Yan,  and  Liu:  Phase-shifts-based  method  of  breast  density  measurement 


4244 


4244 

11  A.  Snigirev  et  al “On  the  possibilities  of  x-ray  phase  contrast  micro¬ 
imaging  by  coherent  high-energy  synchrotron  radiation,”  Rev.  Sci.  Instrum. 
66,5486-5492(1995). 

12 S.  Wilkins,  T.  Gureyev,  D.  Gao,  A.  Pogany,  and  A.  Stevenson,  “Phase  con¬ 
trast  imaging  using  polychromatic  hard  x-ray,”  Nature  (London)  384,  335- 
338  (1996). 

13 F.  Arfelli  et  al.,  “Mammography  with  synchrotron  radiation:  Phase- 
detection  techniques,”  Radiology  215,  286-293  (2000). 

14E.  Donnelly  and  R.  Price,  “Effect  of  kVp  on  edge-enhancement  index  in 
phase-contrast  radiography,”  Med.  Phys.  29,  999-1002  (2002). 

15  X.  Wu  and  H.  Liu,  “Clinical  implementation  of  phase  contrast  x-ray  imag¬ 
ing:  Theoretical  foundation  and  design  considerations,”  Med.  Phys.  30, 
2169-2179  (2003). 

16G.  Hammerstein,  D.  Miller,  D.  White,  M.  Masterson,  H.  Woodard,  and 
J.  Laughlin,  “Absorbed  radiation  dose  in  mammography,”  Radiology  130, 
485—491  (1979). 

17 K.  Nugent,  T.  Gureyev,  D.  Cookson,  D.  Paganin,  and  Z.  Bamea,  “Quan¬ 
titative  phase  imaging  using  hard  x  rays,”  Phys.  Rev.  Lett.  77,  2961-2965 
(1996). 


i8L.  Allen  and  M.  Oxley,  “Phase  retrieval  from  series  of  images  obtained  by 
defocus  variation,”  Opt.  Commun.  199,  65-85  (2001). 

19P.  Cloetens,  R.  Mache,  M.  Schlenker,  and  S.  Lerbs-Mache,  “Quantitative 
phase  tomography  of  Arabidopsis  seeds  reveals  intercellular  void  network,” 
Proc.  Natl.  Acad.  Sci.  U.S.A.  103,  14626-14630  (2006). 

20N.  A.  Dyson,  X-rays  in  Atomic  and  Nuclear  Physics  (Longman,  London, 
1973). 

21 X.  Wu,  H.  Liu,  and  A.  Yan,  “X-ray  phase-attenuation  duality  and  phase 
retrieval,”  Opt.  Lett.  30,  379-381  (2005). 

22  A.  Yan,  X.  Wu,  and  H.  Liu,  “Robustness  of  phase  retrieval  methods  in  x-ray 
phase  contrast  imaging:  A  comparison,”  Med.  Phys.  38,  5073-5080  (201 1). 

23 M.  Bom  and  E.  Wolf,  Principle  of  Optics,  6th  ed.  (Pergamon,  Oxford, 
1980). 

24R.  S.  Saunders  Jr.,  and  E.  Samei,  “A  method  for  modifying  the  image  qual¬ 
ity  parameters  of  digital  radiographic  images,”  Med.  Phys.  30,  3006-3017 
(2003). 

25  G.  Mawdsley,  A.  Tyson,  C.  Peressotti,  R.  Jong,  and  M.  Yaffe,  “Accurate 
estimation  of  compressed  breast  thickness  in  mammography,”  Med.  Phys. 
36,  577-586  (2009). 


Medical  Physics,  Vol.  39,  No.  7,  July  2012 


COL  10(12),  121101(2012) 


CHINESE  OPTICS  LETTERS 


December  10,  2012 


Phase  retrieval  for  hard  X-ray  computed  tomography  of 
samples  with  hybrid  compositions 

Huiqiang  Liu  (MM-&)1,  Yuqi  Ren  (#-i*|-)1,  Han  Guo  (?p  4&01,  Yanling  Xue 

Honglan  Xie  Tiqiao  Xiao  (W#-#)1*,  and  Xizeng  Wu  (^^Jf)2** 

1  Shanghai  Institute  of  Applied  Physics,  Chinese  Academy  of  Sciences,  Shanghai  201204,  China 
2 Department  of  Radiology,  University  of  Alabama  at  Birmingham,  Birmingham,  Alabama  35233,  USA 
*  Corresponding  author:  tqxiao@sinap.ac.cn;  **  corresponding  author:  xwu@uabmc.edu 
Received  May  31,  2012;  accepted  June  20,  2012;  posted  online  September  14,  2012 


X-ray  tomography  of  samples  containing  both  weakly  and  strongly  absorbing  materials  are  necessary  in 
material  and  biomedical  imaging.  Extending  the  validity  of  the  phase-attenuation  duality  (PAD)  method, 
the  propagation-based  phase-contrast  computed  tomography  (PPCT)  of  a  sample  with  hybrid  composi¬ 
tions  of  both  the  light  and  dense  components  with  60  keV  of  synchrotron  radiation  is  investigated.  The 
experimental  results  show  that  the  PAD-based  PPCT  is  effective  in  imaging  both  the  weakly  and  strongly 
absorbing  components  simultaneously.  Compared  with  the  direct  PPCT  technique,  the  PAD-based  PPCT 
technique  demonstrates  its  excellent  capability  in  material  discrimination  and  characterization.  In  addi¬ 
tion,  the  PAD-based  PPCT  exhibits  a  striking  performance  on  the  image  contrast  enhancement  and  noise 
suppression.  Therefore,  this  technique  is  useful  for  material  and  biomedical  imaging  applications,  espe¬ 
cially  when  the  radiation  dose  involved  imposes  a  serious  constraint. 
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X-ray  tomography  of  samples  containing  both  weakly 
and  strongly  absorbing  materials  is  necessary  in  ma¬ 
terial  and  biomedical  imaging!1,2!.  Propagation-based 
phase-contrast  computed  tomography  (PPCT)  is  a  good 
choice  for  this  purpose  because  it  utilizes  highly  sen¬ 
sitive  phase  contrast  and  requires  a  relatively  simple 
setup!3-5!.  Two  different  types  of  approaches  are  avail¬ 
able  for  PPCT.  One  type  is  direct  PPCT,  in  which  the 
tomographic  reconstruction,  such  as  that  based  on  the 
filtered  back-projection  (FBP)  algorithm,  is  directly  ap¬ 
plied  to  the  acquired  angular  projections  with  mixed  at¬ 
tenuation  and  phase  contrast.  The  other  type  is  phase 
retrieval-based  PPCT,  in  which  the  acquired  projections 
are  preprocessed  with  the  phase  retrieval  before  tomo¬ 
graphic  reconstruction.  The  phase  retrieval  on  the  ac¬ 
quired  angular  projections  can  enhance  image  contrast, 
reduce  artifacts,  and  enable  quantitative  characterization 
of  the  sample!6-8!.  However,  multiple  sets  of  projections 
must  be  obtained  at  different  sample-detector  distances 
(SDDs)  to  proceed  with  the  phase  retrieval,  indicating 
that  a  laborious  acquisition  process  is  needed,  and  con¬ 
siderably  more  radiation  dose  is  introduced!9-11!.  There¬ 
fore,  a  phase  retrieval-based  PPCT  method  that  uses 
only  a  single  set  of  projections  is  highly  desirable. 

The  phase  shift  results  from  X-ray  coherent  scattering, 
whereas  the  attenuation  arises  from  the  three  kinds  of 
interactions  between  X-ray  and  matter,  such  as  photo¬ 
electric  absorption,  incoherent  scattering,  and  coherent 
scattering.  The  phase  attenuation  duality  (PAD)  rela¬ 
tion  holds  for  X-rays  of  60-500  keV  transmitted  from  the 
samples  of  light  elements  with  atomic  numbers  Z  <  10!'!. 
In  these  cases,  the  attenuations  and  phase  shifts  of  the 
samples  are  all  related  to  the  projected  electron  densi¬ 
ties  along  the  X-ray  paths.  Under  PAD  conditions,  the 
X-ray  intensity  propagation  equation  is  simplified,  such 
that  robust  phase  retrieval  can  be  achieved  from  only  a 


single  projection  at  each  view  angle.  This  letter  reports 
the  experimental  results  for  a  sample  with  hybrid  com¬ 
position  to  demonstrate  the  capability  of  the  PAD-based 
PPCT  technique  in  material  discrimination  and  charac¬ 
terization. 

A  sample  for  phase-contrast  imaging  is  modeled  as  a 
2D  transmission  function  in  the  object  plane  X=  ( x ,  y), 
with  the  amplitude  denoted  by  A(X)  and  the  phase  by 
(,c(X) .  Maps  A2(X)  and  <^(X)  are  the  attenuation  image 
of  the  sample  and  the  phase  image,  respectively.  A2(X) 
is  related  to  the  projected  electron  density  pe,  P(X)  when 
the  PAD  condition  is  satisfied,  given  that  A2(X)=exp[ 

<7 KNPe,  p(X)],  where  ctkn  is  the  total  cross  section  of 
the  X-ray  Compton  scattering  from  a  single  free  elec¬ 
tron.  However,  the  phase  shift  y>(X)  along  a  ray  path  is 
given  by  Are(oe,  P(X),  and  re  is  the  classic  electron  radius. 
Based  on  the  paraxial  Fresnel-Kirchhoff  diffraction,  the 
projected  electron  density  pe,  P(X)  is  given  as!7! 


Pe,p(X> 


r1 


/{M27(MX;R1+R2)} 
Jin(l  +  2tt(X2^)m2) 


(1) 

where  M  =  (i?i  +  R,2)/Ri,  Ri  is  the  source-to-sample 
distance,  J?2  denotes  the  SDD,  J(X;  R\  +  i?2)  is  the  pro¬ 
jection  intensity  at  detector  plane,  and  is  the  entrance 
intensity.  In  Eq.  (1),  /{  }  denotes  the  2D  Fourier  trans¬ 
form,  /-1{  }  is  its  inverse,  and  u  denotes  the  spatial 
frequency  vector  in  the  object  plane.  The  initial  angular 
projections  in  PAD-based  PPCT  are  preprocessed  with 
PAD-based  phase  retrievals,  and  then  the  standard  FBP 
algorithm  is  employed  to  reconstruct  the  distribution  of 
the  electron  densities  in  the  sample  from  the  prepro¬ 
cessed  projections. 

The  experiments  were  performed  at  the  X-ray  Imaging 
and  Biomedical  Application  Beamline  (BL13W1)  of  the 
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Fig.  1.  Schematic  of  the  propagation-based  imaging  system. 

Shanghai  Synchrotron  Radiation  Facility  (SSRF),  which 
can  provide  high  flux  X-rays  in  the  energy  range  of 
8-72.5  keV.  The  sample  was  located  at  approximately 
34-m  downstream  of  the  wiggler  source  to  provide  both 
high  spatial  coherence  and  almost  parallel  beam.  The 
schematic  of  the  experimental  setup  is  demonstrated  in 
Fig.  1.  The  experiments  were  conducted  at  60  keV, 
with  a  SDD  of  80  cm.  The  phantom  for  the  experi¬ 
ments  was  composed  of  four  types  of  standard  materials 
(GoodFellow  Cambridge  Limited),  namely,  an  alumina 
rod  (AI2O3,  <P  =  1  mm,  3.9  g/cm-3),  an  aluminum  tube 
(Al,  $ln=  1.166  mm,  <£out=  1-42  mm,  2.7  g/cm-3),  a 
PTFE  tube  (C2F4,  <P jn  =  1  mm,  <A,ut=  1-82  mm,  2.2 
g/cm-3),  and  a  PMMA  rod  (C5H8O2,  $  =  1  mm,  1.19 
g/cm-3).  The  metallic  oxide  and  metal  components  were 
chosen  to  simulate  dense  materials,  such  as  the  bones, 
whereas  the  polymer  components  were  used  to  simulate 
light  materials,  such  as  soft  tissues. 

A  total  of  900  phase-contrast  radiographs  were  col¬ 
lected  over  a  180°  rotation  for  a  tomographic  data  set; 
each  radiograph  was  acquired  with  an  exposure  time  of 
208  ms.  The  detector  was  a  CCD  (Photonic  Science) 
with  a  9-/um  effective  pixel  size.  Considering  that  a 
highly  coherent  X-ray  beam  and  large  SSD  were  em¬ 
ployed,  all  these  radiographs  exhibited  not  only  the 
usual  attenuation-based  image  contrast,  but  also  en¬ 
hanced  material  interfaces  as  a  result  of  the  refraction 
and  diffraction  of  the  phase-shifted  X-rays.  For  compar¬ 
ison,  tomograms  were  reconstructed  using  the  angular 
projection  data  set  through  two  different  techniques, 
namely,  the  direct  and  the  PAD-based  PPCTs.  The  to¬ 
mogram  by  direct  PPCT,  in  which  the  standard  FBP  al¬ 
gorithm  was  directly  applied  to  the  initial  projections,  is 
shown  in  Fig.  2(a).  The  structures  for  the  dense  and  light 
material  components  can  all  be  distinguished  from  each 
other  because  of  the  edge  enhancements  at  the  interfaces. 
Notably,  some  micropores  inside  the  AI2O3  rod  became 
apparent  as  the  results  of  the  edge-enhancement.  This 
feature  indicates  that  the  direct  PPCT  method  is  suitable 
for  distinguishing  the  small  cracks  or  voids  embedded  in¬ 
side  a  sample.  However,  the  bulk  contrast  between  the 
different  materials  appeared  rather  weak.  Moreover,  Fig. 
2(a)  presents  anomalously  large  or  even  negative  appar¬ 
ent  linear  attenuation  coefficients  at  material  interfaces, 
impeding  proper  material  characterization.  According 
to  a  general  analysis  of  the  direct  PPCT  technique,  the 
pixel  value  indicated  in  Fig.  2(a),  namely  the  recon¬ 
structed  apparent  attenuation  coefficient,  was  not  equal 
to  the  material  linear  attenuation  coefficient.  Rather, 
the  reconstructed  apparent  attenuation  coefficient  is  a 
sum  of  three  parts,  i.e. ,  the  linear  attenuation  coefficient 
of  the  material,  the  Laplacian  of  the  refractive  index 
of  the  material,  and  the  nonlocal  contribution,  which 
is  dependent  on  the  global  variations  in  the  phantom 


attenuation  coefficients  and  refractive  indices!5/  The 
presence  of  an  excessive  noise  added  uncertainty  to  the 
apparent  attenuation  coefficients. 

The  tomogram  shown  in  Fig.  2(b)  was  reconstructed 
using  the  PAD-based  PPCT  method,  in  which  the  ac¬ 
quired  projections  were  preprocessed  with  the  PAD- 
based  phase  retrievals  to  retrieve  the  pe,  P-maps  for  each 
individual  angular  view.  FBP-based  algorithm  was  then 
applied  to  these  retrieved  pe.  p-maps  for  tomography  re¬ 
construction.  For  comparison,  the  intensity  profiles  along 
two  white  dash-dot  circles  in  Figs.  2(a)  and  (b)  are  shown 
in  Fig.  3.  The  red  dashed  line  associated  with  the  direct 
PPCT  method  shows  severe  high-frequency  noise  and 
weak  bulk  contrast,  which  cannot  provide  useful  means 
for  material  characterization.  However,  the  blue  solid  line 
with  the  PAD-based  PPCT  method  demonstrates  good 
bulk-contrast  and  high  signal-to-noise  ratios,  enabling 
qualitative  and  quantitative  material  characterization. 
The  bulk  contrast-to-noise  ratios  (CNRs)  achieved  in 
Figs.  2(a)  and  (b)  were  compared  to  quantify  the  noise- 
suppression  advantage  of  the  PAD-based  PPCT.  The 
signal  was  defined  as  the  mean  image  intensity  in  a  re¬ 
gion  of  interest,  excluding  the  boundaries.  The  contrast 
was  defined  as  the  difference  in  the  signals  in  two  re¬ 
gions  of  interest,  and  the  noise  was  the  average  of  the 
standard  deviations  of  the  intensities  in  the  two  regions. 
The  bulk  CNRs  achieved  with  the  two  techniques  are 
shown  in  Fig.  4,  indicating  that  the  bulk  CNRs  with  the 
PAD-based  PPCT  method  were  approximately  10  to  15 
times  higher  than  that  with  the  direct  PPCT  method. 
This  significant  enhancement  in  the  bulk  CNR  implies 
that  the  PAD-based  method  has  the  potential  to  reduce 
tremendously  the  radiation  dose.  This  feature  is  partic¬ 
ularly  useful  for  biomedical  applications,  where  a  good 
bulk  contrast  is  critical  for  tissue/lesion  characteriza¬ 
tion.  Some  ring  artifacts  are  present  in  both  Figs.  2(a) 
and  (b),  which  are  due  to  the  beam  quality  degradation 
at  high  photon  energies. 

In  addition,  the  histograms  derived  from  Figs.  2(a) 
and  (b)  are  demonstrated  in  Figs.  5(a)  and  (b),  respec¬ 
tively.  Figure  5(a),  which  is  derived  from  direct  PPCT, 
shows  that  the  overlapped  histogram  peaks  hamper  the 
discrimination  between  the  four  different  materials  in  the 
sample.  However,  Fig.  5(b)  from  the  PAD-based  PPCT 
shows  that  the  well-separated  peaks  clearly  discriminate 
among  the  four  different  materials  in  the  sample.  The 
histogram  peaks  in  Fig.  5(b)  were  identified  as  those 
of  PMMA,  PTFE,  Al,  and  AI2O3  from  left  to  right. 

(cirr1)  x  1023 


Fig.  2.  Reconstructions  for  the  sample  at  60  keV.  Tomograms 
of  (a)  direct  and  (b)  PAD-based  PPCTs. 
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Pixel 

Fig.  3.  (Color  online)  Profiles  of  rescaled  intensity  on  the 
trace  of  white  dash-dot  circle  in  Figs.  2(a)  and  (b),  respec¬ 
tively. 


Fig.  4.  Plots  of  CNR  for  both  direct  and  PAD-based  PPCTs 
as  well  as  the  gain  of  the  CNRs  of  the  PAD-based  and  direct 
PPCTs. 


Fig.  5.  Histograms  of  the  electron  density  distribution  of  the 
samples  with  (a)  direct  and  (b)  PAD-based  PPCTs. 

This  layout  of  the  peaks  is  consistent  with  the  true  elec¬ 
tron  density  values  of  these  materials. 

Moreover,  some  quantitative  information  can  also  be 
obtained  from  the  histogram  shown  in  Fig.  5(b).  The 
PAD  condition  at  60  keV  holds  for  the  materials  made 
of  light  elements  with  Z<  10.  Hence,  good  accuracies 
can  be  achieved  in  recovering  the  electron  densities 
of  the  PTFE  tube  and  PMMA  rod  in  the  sample  be¬ 
cause  these  two  materials  consist  of  the  light  elements 
with  Z<  10.  The  electron  density  of  a  material  is  the¬ 
oretically  expressed  as  pe  =  NapY^[w\{Zi/ Ai)],  where 
Na  is  the  Avogadro  constant  and  p  is  the  mass  den¬ 
sity  of  the  material;  Wi,  Z\,  and  A;  are  the  elemen¬ 
tal  weight  fraction,  atomic  number,  and  atomic  weight 


of  the  material,  respectively.  Therefore,  the  obtained 
theoretical  electron  density  values  were  pe(PMMA)  = 
3.865xl023  cm"3  and  pe(PTFE)  =  6.358xl023  cm"3, 
and  their  ratio  was  pe(PTFE)/pe(PMMA)=1.645.  How¬ 
ever,  the  histogram  in  Fig.  5(b)  shows  that  the  ra¬ 
tio  [/oe(PTFE)/pe(PMMA)]Peak  =  1.714  differs  from  its 
theoretical  value  by  only  4.2%.  Moreover,  Fig.  2(b) 
shows  that  the  reconstructed  electron  density  values  of 
these  two  materials  are  less  than  12%  different  from 
their  theoretical  values.  These  results  demonstrate  that 
the  PAD-based  PPCT  method  achieved  good  accura¬ 
cies  in  evaluating  the  electron  densities  of  the  materi¬ 
als  composed  of  light  elements.  Similar  ratios  involv¬ 
ing  aluminum  atoms,  such  as  [/9e(Al)/pe(PMMA)]peak 
and  [pe(Al203)/pe(PMMA)]peak,  differ  from  their  the¬ 
oretical  ratios  by  as  much  as  36%.  These  larger 
errors  are  expected  considering  that  the  PAD  does 
not  hold  for  heavier  elements,  such  as  aluminum 
(Z=  13).  These  errors  associated  with  heavier  ele¬ 
ments  can  be  reduced  as  the  X-ray  photon  energy 
increases. 

Although  several  other  methods  that  use  only  a  single 
set  of  projections  for  the  phase  retrieval-based  PPCT 
are  available,  these  methods  are  not  applicable  to  sam¬ 
ples  of  hybrid  materials  with  both  light  and  dense 
components.  As  analyzed  in  Ref.  [13],  some  of  these 
methods  require  samples  with  diminished  attenuation, 
whereas  others  are  applicable  only  to  samples  composed 
of  single  or  binary  materials  with  known  linear  attenu¬ 
ation  coefficients  and  refractive  indices.  The  proposed 
method  is  therefore  advantageous  because  it  does  not 
require  such  assumptions  and  constraints.  The  PAD- 
assumption  holds  only  if  the  X-ray  Compton  scatter¬ 
ing  dominates  in  the  X-ray  attenuation  process,  as  we 
point  out  at  the  beginning.  Hence,  additional  research 
is  needed  to  investigate  whether  the  PAD-based  PPCT 
method  is  applicable  for  other  combinations  of  light  and 
dense  materials  encountered  in  material  and  biomedical 
applications. 

In  conclusion,  extending  the  validity  of  the  PAD-based 
method,  the  PPCT  of  a  sample  with  hybrid  composi¬ 
tions  is  investigated  using  a  60-keV  synchrotron  radia¬ 
tion.  The  experimental  results  show  that  the  PAD-based 
PPCT  is  effective  in  imaging  simultaneously  both  the 
weakly  and  strongly  absorbing  components.  Compared 
with  the  direct  PPCT  technique,  the  PAD-based  PPCT 
demonstrats  its  excellent  capability  in  material  discrim¬ 
ination  and  characterization.  This  PAD-based  method 
also  exhibits  a  striking  performance  on  the  bulk-contrast 
enhancement  and  noise  suppression.  This  technique 
therefore  has  great  potential  for  material  and  biomedical 
applications,  especially  when  the  radiation  dose  involved 
imposes  a  serious  constraint. 
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